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C  dimensionless  number,  C  / C 

sw  ’  s'  w 

cy  specific  heat  at  constant  volume,  J/kg-K 

D  diagonal  distance,  m 

E  enthalpy,  J/kg 

F  flow  rate  through  a  control- volume  face,  kg/s 
H  latent  heat,  J/kg,  Sections  IV  and  V 
H  heat  transfer  coefficient,  V/(m  -K),  Section  III 
h  enthalpy,  J/kg,  Sections  I  and  II 

h  height  of  the  spot  heated  heat  pipe,  m,  Section  III 

h  total  thickness,  6  +  6,  m.  Section  V 

’  p  7  ’ 

* 

h  dimensionless  number,  h/R^ 


xiv 


Nomenclature  (continued) 


hr 

fg 

h- 

1 

K 

k 

k£w 

ksw 

Kvl 

Klv 

L 

L 

LC 

LE 

lh 

L1,L2 

i 

m 

M 

\ 

N 

r 

■t 

N 


latent  heat  of  evaporation,  J/kg 

o 

heat  transfer  coefficient  at  interface,  V/m  -K 

ratio  of  the  thermal  conductivities  of  the  wall  to  the  fluid, 

Vkf 

thermal  conductivity,  V/m-K 

dimensionless  number,  k^/ky 

dimensionless  number,  k  /k 

s  w 

solid  wall  and  effective  liquid  thermal  conductivity  ratio, 

Vkeff 

effective  liquid  and  vapor  thermal  conductivity  ratio, 

Wkv 

total  length  of  the  pipe,  m,  Sections  I  and  II 

reference  length,  m,  Section  IV 

condenser  length,  m 

evaporator  length,  m 

length  of  heated  block,  m 

length  of  spot  heated  heat  pipe,  m 

interface  position  along  the  diagonal,  m,  Sections  IV  and  V 
mass  flux,  kg/m  - s 
Mach  number 

dimensionless  number,  q^R^/ (Tm_ T^)kw 

9 

dimensionless  number,  R,  (T  -T  )  at/ k 

’  h v  m  i'  '  w 

dimensionless  number,  T - / (T  -T) 

’  i' '  m 

dimensionless  number,  UC  (T  -T)k 

’  wv  m  i;  w 
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Nomenclature  (continued) 

Nuz  local  Nusselt  number  at  the  solid- fluid  interface,  2r^h^/k| 

n  normal  outward  direction 

* 

n  dimensionless  normal  outward  direction,  n/R^ 

o 

P  pressure  of  fluid,  N/m 

Pq  datum  pressure  at  the  end  cap  of  the  evaporator  at  the 

2 

liquid- vapor  interface,  N/m 
Pr  Prandtl  number,  pCp/k 

Pe  Peclet  number,  PrRe 

Q  heat  source  power,  V,  Sections  IV  and  V 

Q  total  heat  input  rate  at  evaporator  outer  pipe  wall,  V 

Section  I 

Qg  input  heat,  V 

radiation  heat  transfer  to  surroundings,  V 
2 

q  heat  flux,  M/m 

2 

qg  heat  flux,  V/m 

R  radius  of  heat  pipe,  m,  Section  III 
R  gas  constant  in  Eq.(1.5a),  J/kg-K 

Re  inlet  axial  Reynolds  number,  2r^ 

Ren  axial  Reynolds  number  at  the  adiabatic  section,  2R  v  / v 
a  v  a' 

Rei  radial  Reynolds  number  at  the  solid- fluid  interface,  v-r./i/ 
R^  heat  source  radius,  m 

Ro  radius  of  the  numerical  domain,  m 

Ry  vapor  space  radius,  m 
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Nomenclature  (continued) 

outer  pipe  wall  radius,  m 

Rer  average  radial  Reynolds  number  at  the  liquid- vapor 
interface,  p  v  R  la 
r  radial  coordinate,  m 

S  coefficient  in  Eq.  (4.7)  or  Eq.  (5.19) 
source  term  in  the  energy  equation 
St  Stefan  number,  c  (T  -T  )/H  or  c  (T  -T.)/H 
s  circumferential  distance  coordinate,  m,  Section  III 

s  dimensionless  interface  position  along  the  diagonal,  £/ D 
Section  IV 

SH  Vjj/2tR 

T  temperature,  K 

T*  (T  -  Ts)/Th,  Section  III 

* 

T  dimensionless  temperature,  (T- T^)/(T  -  T^) ,  Section  V 

T+  The  "Kirchoff"  temperature 

Tjj  reference  temperature,  K 

initial  temperature,  K 

Tm  melting  temperature,  K,  Sections  IV  and  V 
Tq  temperature  of  environment,  K 

TQ  temperature  of  evaporator  end  cap  at  liquid- vapor  interfa 

K 

T  vapor  temperature,  K 

u 

wall  temperature,  K 
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Nomenclature  (continued) 

the  temperature  defined  in  Fig.  4.2,  K 
T2  the  temperature  defined  in  Fig.  4.2,  K 
t  time,  s 

t  dimensionless  number,  tkw/(CwR^  ) 

AT  T  -  T  ,  °C,  Section  III 

AT  temperature  range  T2  -  T^,  K,  Sections  IV  and  V 

ATc  critical  temperature  difference,  K 

At  time  step,  s 

U  heat  source  velocity,  m/s 

* 

U  dimensionless  velocity,  C^UR^/k^ 

u,v,w  velocities,  m/s 

Vq  evaporation  or  condensation  velocity  at  the  vapor- liquid 
interface,  m/s 

vn  normal  melting/freezing  speed,  m/s 

3 

AV  volume  of  the  control  volume,  m 

wa  average  axial  velocity  in  the  adiabatic  section,  m/s 

Vg  width  of  heated  block,  m 

wz  average  local  axial  velocity  over  the  cross  section  at  any 
axial  location,  m/s 
width  of  spot  heated  pipe,  m 
x  x/tR,  Section  III 

Xj  distance  of  the  edge  of  heated  block  from  the  origin,  m 
x,y,z  moving  coordinate  directions  in  Section  V  and  general  coordinate 
directions  in  other  sections 
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Nomenclature  (continued) 

x',y',z'  fixed  coordinate  directions,  Section  IV 
*  *  * 

x  ,  y  ,  z  dimensionless  distances,  x/R^,  y »  z/\ 

X-X  the  plane  indicated  in  Fig.  5.4  and  Fig.  5.14 

Y  s/tR 

Greek  Symbols 

2 

a  thermal  diffusivity,  m  /s 

F  coefficient  in  Eq.  (4.7)  or  (5.19)  and  exchange  coefficient  k, 
Sections  I  and  II 

7  ratio  of  specific  heats,  cp/cy 
b  wall  thickness,  m,  Sections  III,  IV  and  V 

5  wall  or  liquid- wick  thickness,  m,  Section  I 
* 

6  dimensionless  wall  thickness,  <5/R^ 

<5m  maximum  melting  front  depth,  m 
* 

<5m  dimensionless  maximum  melting  front  depth,  <5m/R^ 

<5p  thickness  of  the  PCM,  m 
* 

b  dimensionless  thickness  of  the  PCM,  b  /R, 

p  ’  p'  h 

bx  x- direction  distance  between  two  adjacent  grid  points 

A  dimensionless  wall  or  liquid-wick  thickness,  b/ 2Ry,  Section 

I;  (rQ  -  r-)/2r-,  Section  II 

Ax  x- direction  width  of  the  control  volume 

by, by  similar  to  Ax,  bx 

Az,<5z  similar  to  Ax,  bx 
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Nomenclature  (continued) 
r]  dummy  variable 

f  emissivity 

e  wick  porosity,  Sections  I  and  II 

0  dimensionless  temperature,  Section  I  and  II;  coordinate 

direction,  Sections  IV  and  V 

0 ^  dimensionless  temperature,  (T  -Tm)/(Tm  -T^) 

0Q  dimensionless  temperature,  (T^  -  T  ) / (T  -  T  ) 

p  dynamic  viscosity  of  fluid,  kg/m-  s 

3 

p  density  of  fluid,  kg/m 

<ji  viscous  dissipation  term  in  the  energy  eqn.  (1.4),  Section  I; 

dummy  variable,  Section  IV 

2 

r  dimensionless  time,  ta/L 

2  4 

a  Stefan- Boltzmann  constant,  V/(m  -K  ) 

Subscripts 
a  adiabatic 

B  "Bottom"  neighbor  of  grid  P 

b  control- volume  face  between  P  and  B,  Sections  IV  and  V 
b  bulk  condition,  Sections  I  and  II 

C  condensation 

c  condenser 

E  evaporation 

E  "East"  neighbor  of  grid  P,  Sections  IV  and  V 
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Nomenclature  (continued) 

e  control  volume  face  between  P  and  E,  Sections  IV  and  V 
e  evaporator 

eff  effective 

f  fluid 

h  harmonic  mean 

i  initial  condition,  Sections  IV  and  V 

i  vapor- liquid  interface  for  Section  I,  solid  fluid  interface  for 

Section  I 

in  entrance  condition  to  the  pipe 
1  liquid 

lw  liquid- wick 

m  mean  value  in  Section  I  and  II 
m  mushy  pahse  in  Section  IV  and  V 
N  "North"  neighbor  of  grid  P 
n  contol- volume  face  between  P  and  N 

o  variable  at  outer  pipe  wall 

P  grid  point 

r  radial  direction 

S  "South"  neighbor  of  grid  P 

s  control- volume  face  between  P  and  S,  Sections  IV  and  V 
T  "Top"  neighbor  of  grid  P 
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Nomenclature  (continued) 
t  control- volume  face  between  P  and  T 

v  vapor 

V  "Vest"  neighbor  of  grid  P 
w  control- volume  face  between  P  and  V 

w  wall 

z  at  z  location 

0  variable  with  initial  value  before  iteration  and  the 

properties  corresponding  to  Tq 

Superscripts 

+  dimensionless  variable 

average  variable 
*  dimensionless  variable 
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I.  NUMERICAL  ANALYSIS  OF  THE  VAPOR  FLOV  AND  THE  HEAT  CONDUCTION 
THROUGH  THE  LIQUID- VICK  AND  THE  WALL  IN  A  HEAT  PIPE 

1 . 1  Summary 

A  numerical  analysis  is  presented  for  the  overall  performance  of  heat 
pipes  with  single  or  multiple  heat  sources.  The  analysis  includes  the  heat 
conduction  in  the  wall  and  liquid- wick  regions  as  well  as  the 
compressibility  effect  of  the  vapor  flow  inside  the  heat  pipe.  The 

two-dimensional  elliptic  governing  equations  in  conjunction  with  the 
thermodynamic  equilibrium  relation  and  appropriate  boundary  conditions  are 
solved  numerically  for  the  whole  domain.  The  solutions  are  compared  with 
existing  experimental  data  for  the  vapor  and  wall  temperatures  at  both  low 
and  high  operating  temperatures.  The  effects  of  conjugate  axial 
conduction,  vapor  compressibility,  flow  reversal,  and  viscous  dissipation 
are  discussed  as  well  as  the  accuracy  of  the  parabolic  version  versus  the 
elliptic  version.  The  results  show  that  the  axial  wall  conduction  tends  to 
distribute  the  temperature  more  uniformly  for  the  heat  pipe  with  large 
solid  wall  and  effective  liquid  thermal  conductivity  ratios.  The 

compressible  and  incompressible  models  show  a  very  close  agreement  for  the 
total  pressure  drop  while  the  local  pressure  variations  along  the  heat  pipe 
are  quite  different  for  these  two  models  for  the  cases  with  high  radial 
Reynolds  numbers  at  the  interface.  Also,  the  partially  parabolic  solution 
provides  fairly  accurate  results  compared  with  the  elliptic  solution  except 
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for  the  cases  with  a  large  radial  Reynolds  number  at  the  interface  or  large 
solid  wall  and  effective  liquid  thermal  conductivity  ratios. 
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1.2  Introduction 


The  heat  pipe  is  a  device  which  has  a  very  high  thermal  conductance. 
It  is  a  closed  vacuum  tube  or  chamber  of  different  shapes  whose  inner 
surfaces  are  lined  with  a  porous  capillary  wick  (Fig.  1.1).  The  wick  is 
saturated  with  the  liquid  phase  of  a  working  fluid,  and  the  remaining 
volume  of  the  tube  contains  the  vapor  phase.  Heat  applied  at  the 
evaporator  causes  the  liquid  to  vaporize  and  the  vapor  to  move  to  the  cold 
end  of  the  tube  where  it  is  condensed.  The  condensate  is  returned  to  the 
evaporator  by  capillary  forces.  The  amount  of  heat  that  can  be  transported 
as  the  latent  heat  of  vaporization  is  usually  several  orders  of  magnitude 
larger  than  that  which  can  be  transported  in  a  conventional  convective  heat 
transfer  system.  Thus,  a  small  heat  pipe  can  transport  a  large  amount  of 
heat. 

The  first  recorded  patent  of  a  device  similar  to  the  heat  pipe  was  by 
Perkins  and  Buck  (1892).  It  is  called  a  Perkins  tube  which  is  actually  a 
wickless  heat  pipe.  It  relies  on  gravitational  forces  to  return  the  liquid 
from  the  condenser  to  the  evaporator.  In  1942,  Gaugler  applied  a  wick 
structure  to  the  Perkins  tube  and  invented  the  heat  pipe  which  can  operate 
at  different  orientations.  However,  it  was  not  widely  publicized  until 
1964  when  Grover  et  al.  at  the  Los  Alamos  Scientific  laboratory 
independently  reinvented  the  concept.  Grover  (1963)  included  a  theoretical 
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Figure  1,1  The  Multiple  Evaporator  Heat  Pipe 
and  Coordinate  System, 


analysis  and  experimental  results  in  his  patent.  The  demonstration  heat 
pipe  was  made  from  a  stainless  steel  pipe  incorporating  a  wire  mesh  wick 
and  sodium  as  the  working  fluid,  although  lithium  and  silver  were  also 
mentioned  as  possible  working  fluids.  Since  then,  the  heat  pipe  has  gained 
ever  increasing  development  and  applications. 

The  performance  of  heat  pipes  is  restricted  by  several  limitations, 
namely,  (1)  vapor  flow  limit:  sonic  and  viscous  limits;  (2)  liquid  flow 
limit:  capillary,  entrainment,  and  boiling  limits.  The  sonic  limit  is 
defined  as  when  the  vapor  velocity  is  so  high  that  it  reaches  the  speed  of 
sound  (»=i),  thus  affecting  the  overall  heat  pipe  performance  because  of 
poor  temperature  uniformity.  The  viscous  limit  usually  occurs  at  low- 
working  temperatures  where  the  viscous  forces  are  dominant  in  the  vapor 
flow.  As  described  by  Busse  (1973),  the  maximum  heat  flux  due  to  the 
viscous  limit  occurs  when  the  vapor  pressure  is  reduced  to  zero  (P=0).  The 
capillary  limit  is  the  wicking  limitation  whenever  the  capillary  forces  are 
insufficient  to  overcome  the  vapor  and  liquid  pressure  losses  (e.g.,  APc  < 
AP^  +  A?v ,  where  APC  is  the  capillary  pressure,  and  AP^  and  APy  are  the 
liquid  and  vapor  pressure  losses,  respectively).  The  entrainment  limit 
represents  an  interaction  between  the  countercurrent  vapor  and  liquid  that 
leads  to  the  entrainment  of  liquid  droplets  from  the  wick  structure,  which 
leads  to  an  insufficient  amount  of  liquid  returning  to  the  evaporator  which 
limits  the  heat  pipe  operation.  The  boiling  limit  is  a  limitation  of  the 
radial  heat  flux  density.  For  a  considerably  high  ladial  heat  flux  at  the 
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evaporator,  vapor  bubbles  may  be  formed  in  the  evaporator  wick.  The 
formation  of  vapor  bubbles  in  the  wick  structure  can  cause  hot  spots  and 
obstruct  the  circulation  of  the  liquid. 

A  common  type  of  heat  pipe  which  meets  the  growing  needs  in  space 
applications  (e.g.,  transportation  of  energy  in  the  reactor  core  of  nuclear 
power  systems  (Koening,  1985)  and  supersonic  jet  surface  cooling 
(Silverstein,  1971))  is  a  liquid  metal  heat  pipe,  which  uses  a  liquid  metal 
for  the  working  fluid  such  as  sodium.  The  major  concern  for  this  type  of 
heat  pipe  is  the  sonic  limit  during  the  starting  period  of  operation,  i.e., 
at  the  lower  working  temperature  range  for  a  certain  working  fluid.  This 
is  because  the  vapor  has  a  low  density  and  a  high  viscosity.  In  this  case, 
the  compressibility  of  the  vapor  and  the  viscous  dissipation  may  have  a 
significant  effect  on  the  performance  of  the  heat  pipe.  Also,  the  high 
vapor  velocity  and  viscosity  may  cause  a  significant  pressure  drop  which  is 
related  to  the  temperature  profile  by  the  thermodynamic  equilibrium  at  the 
vapor- liquid  interface,  thus  affecting  the  overall  heat  pipe  performance. 
The  purpose  of  the  present  study  is  to  determine  the  effects  of  conjugate 
heat  transfer  and  vapor  compressibility  for  liquid  metal  heat  pipes  by 
solving  the  complete  conservation  of  mass,  Navier- Stokes,  and  conservation 
of  energy  equations  by  using  a  computer  code  called  PHOENICS.  The  present 
analysis  includes  the  conjugate  heat  conduction  through  the  wall  and  the 
liquid- wick  matrix,  which  is  usually  neglected  but  is  very  important  for 
heat  pipes  with  multiple  heat  sources.  The  investigation  is  also  extended 
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to  lower  temperature  heat  pipes  such  as  those  using  water  as  the  working 
fluid. 

Literature  Review 


Since  1964,  many  experimental,  analytical,  and  numerical  investigations 
were  made  on  heat  pipes.  The  related  work  presented  here  is  grouped  into 
the  following  three  categories:  (1)  vapor  dynamics,  (2)  coupled  temperature 
field  and  vapor  dynamics,  and  (3)  experimental  studies.  A  brief 
description  of  the  previous  theoretical  and  experimental  work  is  also 
summarized  in  Tables  1.1  and  1.2. 

Vapor  Dynamics 

The  prediction  of  the  vapor  pressure  distribution  in  heat  pipes  was 
first  studied  by  Cotter  (1965)  for  incompressible  flow  and  the  limiting 
cases  of  Rer  <<  1  and  Rer  >>  1.  He  concluded  that  for  Rer  <<  1,  the 
viscous  effects  dominate,  and  the  pressure  gradient  is  the  same  as  in  a 
Poiseuille  flow  while  for  Rer  >>  1,  the  flow  is  dominated  by  inertial 
effects.  The  pressure  recovery  in  the  condenser  was  4/t2  of  the  pressure 
drop  at  the  evaporator. 

Busse  (1967)  further  analyzed  the  hydrodynamics  of  laminar 
incompressible  flow  in  long  cylindrical  heat  pipes.  An  analytical  approach 


7 


«  JZ  -V 
M  V  C 
-  L  -  0 
C  CL  2 

s  *»: 


**  0  0 

•S6-8 

^  U  L 

-•  a  0 

w  a.  o. 


o*o  +■ 
-  s  3 
>  s  «•* 

«  IS 


O  U  O  0 

i  o  i  5. 

—  M  0 


x  v  a 

~  X* 

a  — 

C  0  0 


6  S  “S  t 

<•«  u 

X  4>  u  *U 

at  *j  g*  o  ■-< 

w  —  DL.  3 

—  u  c  cr 
3  o  *  *o  ■•+ 

05  ^  W  4)  *-» 

U  0  4-*  «J 
U  >  C  0 

41  L  « 

—  *  a  «  3  S 

a  v  41  41  ^  & 

g  l  -c  l  e  •« 

—  3  a  w  a 

u  «  g  D 

34)  (a  •-  0  0.  *J 

g  c  g  a  « 

u  0  u  3  g 

z  an  j  *-  -c 


—  y  g  « 

a  o  >«<*-.  • 

C  ri  *— «  y  0 

o  o  •«  ■©  g 

—  u  •  4  i 

—  c  x  c  u 

s  s  g  ~  —  g 

3  <**  l  a  «  n 

tJO  0  X 

CL.*-  •  V  0 


^  h.  *—  K 

T  j{  «  ) 
c  3  #  5 


♦  9  fcS 

*4  W-H 

«  S  Jj  u 

iia 


o 

^  0  •* 
eu*-  — 
u  -  5- 
•h  *j  i  e 

j*;  ?  0 
Is;* 


•5  S  *■? 

o  -  a  n 

—  *»  g  x 

- 

*  =>  -  ® 
a  a  c 
s 

n  .  — 

-  0  T> 

8  9  g  c 
5.x:  o 

;iv, 
o  a>  o  d 
u  >  x 
a.  v  0  n 
a.  -o  44  a 

0  •*»  g 

a  a  0 

C  0  -*4  — 

<*'-*© 


c  ^  c 
0  0  0  0 
3  -  C 

-  0  44  0 

o  g  g  -> 

3  C  C  w 
O*  0  3  0 
0>  *J  W  s 


o  >  o  w 
•--<X  L 
*4  o  **  o 
a  «  > 

2  3  th 

U  )  I  8 


41  II  0  U 

u  5  e  u 

cl  g  0  cl 


l  1 


8 


Popov  (1976)  No  Well  Momentum  (equilibrium  two-  the  eonic  limit  in  terms  of 

phase)  temperature  Mere  presented 

for  high  temperature  liquid 
metal  heat  pipes 
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_ TABLE  1  -P  Snate  Related  Bxperiaental  Work  on  Heat  Pipes _ 
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Faghn  annular  water  annulus  axial  grooves  outer  wall  steady  state  capillary  in 
and  Thomas  heat  pipe  t  emper  at  ore  was  measured 
(1988) 


was  used  to  solve  the  momentum  and  continuity  equations.  The  heat  pipe 
under  consideration  had  three  distinct  sections  with  constant  mass 
injection  in  the  evaporator  and  extraction  in  the  condenser.  The  length  of 
each  section  was  assumed  to  be  much  larger  than  the  diameter  of  the  vapor 
channel.  To  solve  the  momentum  equation,  the  axial  velocity  profile  was 
approximated  by  a  fourth- order  polynomial  with  respect  to  the  radius  and  a 
correction  function.  He  found  that  the  normalized  axial  velocity  component 
is  constant  along  the  evaporator,  approaches  the  Poiseuille  profile  in  the 
adiabatic  section,  and  varies  strongly  along  the  condenser.  The  expression 
for  the  pressure  distributions  along  the  evaporator,  adiabatic,  and 
condenser  sections  were  derived  separately. 

Bankston  and  Smith  (1973)  analyzed  the  incompressible  two-dimensional 
laminar  vapor  flow  in  a  heat  pipe.  The  solution  was  presented  for  the 
pressure  variation  along  the  pipe  for  various  evaporator  and  condenser 
lengths  at  steady- state  operation.  The  finite  difference  scheme  was 
employed  to  solve  the  momentum  equation  in  terms  of  the  stream  function  and 
vorticity.  In  addition,  a  new  series  solution  for  the  slow-motion  case  was 
also  obtained,  and  it  confirms  numerical  results  in  the  limit  of  low 
Reynolds  numbers.  Various  cases  with  different  radial  Reynolds  number  were 
analyzed,  and  they  concluded  that  flow  reversal  in  the  condenser  may  occur 
when  the  condenser  radial  Reynolds  number  exceeds  2. 

Kadaner  and  Rassadkin  (1975)  studied  the  laminar  vapor  flow  along  the 
heat  pipe.  The  problem  was  solved  by  using  a  parametric  method  which 
introduces  a  velocity  profile  in  terms  of  two  unknown  functions.  The 
integrated  momentum  equation  was  then  solved  in  terms  of  these  functions. 
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The  functions  were  determined  by  analyzing  the  limiting  situation,  and  the 
results  were  presented  in  a  graphic  manner  for  the  determination  of  the 
vapor  pressure  loss  along  the  pipe.  The  dimensionless  pressure  drop  and 
friction  coefficient  were  presented  in  terms  of  the  radial  Reynolds  numbers 
and  the  dimensionless  axial  distance. 

Bystrov  and  Mikhailov  (1982)  also  studied  the  laminar  vapor  flow  in  the 
heat  pipe  condenser  using  the  parametric  method.  The  axial  velocity 
profile  was  approximated  by  a  fourth- order  polynomial  with  respect  to  the 
radius  and  an  unknown  parameter  which  was  later  determined  by  the  boundary 
conditions.  The  results  were  presented  for  pressure  and  velocity 
variations  with  the  Reynolds  number  in  the  range  of  1  to  30.  The 
comparison  of  the  theoretical  results  with  the  experimental  data  by  Quaile 
and  Levy  (1972)  showed  a  good  agreement  under  the  conditions  studied. 

Faghri  (1986)  solved  the  2-D  incompressible  laminar  vapor  flow  in  the 
concentric  annular  heat  pipe.  The  parabolic  Navier- Stokes  equation  plus 
the  continuity  equation  were  solved  by  using  the  finite  difference 
numerical  scheme.  The  pressure  drops  and  velocity  profiles  were  presented 
for  symmetric  and  asymmetric  heating  and  cooling  cases.  He  concluded  that 
for  radial  Reynolds  numbers  much  less  than  1,  viscous  effects  dominate,  and 
the  axial  velocity  profile  throughout  the  heat  pipe  is  close  to  the 
Poiseuille  flow  conditions.  For  radial  Reynolds  numbers  greater  than  1, 
the  evaporation  and  condensation  cases  become  qualitatively  different. 
Flow  reversal  was  also  noticed  in  the  condenser  section  for  high 
condensation  cooling  rates.  The  observations  obtained  here  for  an  annular 
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heat  pipe  are  similar  to  what  Busse  (1967)  concluded  for  a  cylindrical  heat 
pipe. 

Narayana  (1986)  studied  the  2-D  steady,  incompressible  laminar  vapor 
flow  in  a  heat  pipe.  The  continuity  and  Navier- Stokes  equations  with  the 
boundary  layer  assumptions  were  solved.  The  effect  of  frictional  and 
inertial  pressure  drops  on  the  static  pressure  is  compared  for  the  values 
of  the  wall  Reynolds  number  in  the  range  of  2  to  -5. 

Busse  (1987a)  presented  a  logical  criterion  for  axisymmetric  internal 
flows  which  indicates  whether  or  not  a  given  velocity  profile  tends  to 
develop  into  a  separation  profile.  He  found  that  the  reversal  of  the 
profile  depends  on  the  profile  shape  and  the  Mach  number.  The  extraction 
of  mass  usually  leads  to  a  much  stronger  profile  variation  than  injection. 
The  viscous  forces  and  mass  injection  tend  to  establish  an  equilibrium 
profile  while  mass  extraction  is  an  inherently  unstable  process  which 
promotes  a  run- away  from  the  equilibrium  condition.  It  is  still  not 
certain,  however,  whether  viscous  forces  can  also  cause  a  flow  reversal. 

Busse  (1987b)  studied  laminar  subsonic  flow  in  cylindrical  heat  pipes. 

The  vapor  was  modeled  as  a  2-D  isothermal  perfect  gas  with  a  constant  heat 

of  vaporization.  The  approach  involves  the  use  of  the  boundary  layer 

approximation  and  a  noncontinuous  power  series  to  describe  the  velocity 

profile.  Based  on  the  analytical  method,  the  numerical  results  for  the 

pressure  recovery  in  the  condenser  were  compared  to  experimental 

measurements.  He  found  that  significant  flow  reversal  at  the  condenser 

wall  first  appears  at  -6.0  >  Re  >  -10.0  giving  rise  to  pressure  recovery 

r .  c 
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in  the  condenser.  Also  confirmed,  is  that  the  onset  of  turbulence  is  at 

Re  «  -6  at  axial  Reynolds  numbers  of  e^y  a  few  hundred.  This 
r,c 

phenomenon  coincides  with  the  appearance  of  significant  reversal  flow.  The 
restriction  for  the  boundary  layer  approximation  is  that  Lc/D  should  not  be 
too  small.  The  results  presented  in  this  report  give  a  fairly  good 

comparison  with  experimental  measurements  with  Lc/D  =4.4.  This  work  is  an 
extension  of  the  work  by  Busse  and  Prenger  (1984)  using  a  similar  model, 
but  the  more  recent  analysis  emphasizes  identifying  the  influence  of 
transition  to  turbulence  on  the  condenser  pressure  recovery. 

Recently,  Bowman  and  Hitchcock  (1988)  modeled  the  unsteady  vapor  flow 
in  a  heat  pipe.  The  vapor  was  assumed  to  be  compressible,  viscous  and 
changing  from  laminar  to  turbulent.  The  Navier- Stokes  equations  were 
solved  numerically  using  the  MacCormack  explicit  finite  difference  method. 
Both  transient  and  steady- state  results  were  obtained.  The  experimental 
data  for  the  pressure  variation  was  obtained  by  using  a  porous  pipe  with 
air  blowing  along  half  of  the  pipe  wall  and  suction  along  the  other  half  of 
the  pipe  which  simulates  the  heat  pipe  vapor  flow.  The  pressure  variation 
of  both  subsonic  and  supersonic  vapor  flow  was  measured  and  compared  with 
the  numerical  results  and  a  fairly  good  agreement  was  found.  The 
experiment  also  demonstrated  that  the  vapor  flow  in  the  heat  pipe  can  be 
accurately  modeled  using  the  steady-  state  governing  equations  because  the 
time  period  of  the  vapor  transient  is  very  short  when  compared  to  the  heat 
pipe  wall  and  wick  transients. 

Faghri  and  Parvani  (1988)  further  analyzed  the  gas  dynamics  of  annular 
heat  pipes.  The  vapor  was  modeled  as  laminar  incompressible  flow.  The 
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elliptic  and  partially  parabolic  Navier- Stokes  equations  were  solved 
numerically  and  compared  for  the  pressure  drop  and  velocity  profile  along 
the  pipe.  For  the  cases  presented,  only  a  small  deviation  was  found 
between  these  two  solutions.  Flow  reversal  was  also  noticed  at  high  radial 
Reynolds  numbers  in  the  condenser  and  long  condenser  lengths.  Based  on  the 
observation  of  a  very  short  hydrodynamic  entry  length  for  the  evaporator 
section,  a  similarity  solution  was  also  proposed  for  the  prediction  of  the 
pressure  drop  and  velocity  profile  for  both  conventional  and  annular  heat 
pipes. 

The  investigators  mentioned  above  only  studied  the  vapor  gas- dynamics 
by  solving  the  mass  and  momentum  conservation  equations.  The  analyses  with 
the  coupling  of  the  energy  equati  will  be  presented  below. 


Coupled  Temperature  Field  and  Vapor  Dynamics 


Levy  (1968)  analyzed  the  one- dimensional  compressible  vapor  flow  in  the 
evaporator  section  of  a  heat  pipe.  He  assumed  a  uniform  temperature 
distribution  and  averaged  the  vapor  velocity  along  tne  radial  direction. 
The  analysis  was  started  from  the  mass,  momentum  and  energy  differential 
equations  incorporating  the  Clausius- Clapeyron  equation,  the  two- phase 
equilibrium  relations  and  the  Clapeyron- Mendeleev  equation  (v  t  =  p--)  for 
the  specific  volume  of  the  saturated  vapor.  The  vapor  quality,  velocity, 
and  pressure  differential  equations  were  obtained  with  respect  to  the  axial 
distance.  These  equations  were  integrated  by  using  the  Runge-Kutta  method, 
and  the  results  under  the  choking  condition  were  presented  for  sodium  heat 
pipes.  Comparing  the  theoretical  calculations  with  existing  data  suggests 
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that  at  sufficiently  low  vapor  pressures,  it  is  the  choking  phenomenon  of 
the  vapor  flow  rather  than  the  capillary  limit  of  the  wick  which  will  limit 
the  maximum  heat  flux  of  the  device.  Bystrov  and  Popov  (1976)  performed  a 
similar  analysis  by  including  the  effect  of  friction  along  the  pipe. 

Rohani  and  Tien  (1973)  studied  the  heat  and  mass  transfer  in  the 
vapor- gas  region  of  a  gas- loaded  heat  pipe.  The  elliptic  governing 
equations  were  solved  in  conjunction  with  the  overall  energy  and  mass 
conservation  balances  and  the  thermodynamic  equilibrium  condition  for  three 
cases  with  different  working  fluids.  They  found  that  in  liquid  metal 
gas- loaded  heat  pipes,  the  vapor- gas  two-dimensional  diffusion  must  be 
considered  in  the  analysis.  A  simple  numerical  framework  to  include  the 
axial  conduction  in  the  wall  and  liquid- wick  matrix  was  also  proposed  in 
that  paper  but  no  results  were  presented. 

In  1974,  Tien  and  Rohani  solved  th  e2- D  elliptic  coupled  momentum  and 

energy  equations  with  the  thermodynamic  equilibrium  relation  at  the 

interface  for  the  heat  pipe  vapor  flow.  The  model  was  laminar, 

compressible  and  steady  state  with  sodium  as  the  working  fluid.  They  used 

a  vorticity- stream- function  approach  similar  to  the  incompressible  flow 

model  used  by  Bankston  and  Smith  (1973).  For  the  boundary  conditions,  they 

specified  that  the  evaporator  and  condenser  are  surrounded  by  corstant  but 

different  ambient  temperatures  and  have  uniform  overall  heat  transfer 

coefficients  between  the  vapor- liquid  interface  and  the  ambient.  The 

pressure  and  temperature  variations  were  presented  in  terms  of  five 

different  radial  Reynolds  numbers:  j Re  |  =  2,  4,  8,  24,  36  or  | Re  |  = 

I*  y  G  r  j  c 

1.33,  2.66,  5.33,  16.0,  24.0.  A  large  deviation  was  found  for  the  elliptic 
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and  parabolic  pressure  drop  results  at  high  radial  Reynolds  numbers.  No 
comparison  was  made  with  experimental  data. 

Bianchi  (1987)  performed  a  similar  analysis  by  using  the  finite  element 
method,  but  he  modeled  the  vapor  flow  as  2-D  and  incompressible  with  water 
as  the  working  fluid  for  the  presentation  of  the  results.  No  pressure 
recovery  in  the  condenser  was  observed  for  the  cases  presented. 

Ismail  et  al.  (1987)  investigated  the  vapor  and  liquid  flows  inside  the 
heat  pipe.  The  finite  difference  method  was  used  to  solve  this  problem. 
The  liquid  and  vapor  were  both  assumed  to  be  2-D  laminar  and 
incompressible.  The  thermodynamic  equilibrium  relation  linked  the 
temperature  and  pressure  at  the  liquid- vapor  interface.  The  liquid  flow  in 
the  porous  medium  is  assumed  to  be  one  phase  and  covers  the  whole  wick 
matrix,  so  the  liquid  flow  was  treated  like  pipe  flow  in  a  porous  annulus 
with  blowing  and  suction  at  the  boundary.  The  overall  pressure  effect  was 
considered  rather  than  considering  the  pressure  difference  at  each 
meniscus.  Flow  reversal  of  the  vapor  was  observed  in  the  evaporator, 
adiabatic  and  condenser  regions,  which  usually  only  happens  in  the 
condenser  region. 

Colwell  et  al.  (1987)  modeled  a  liquid  metal  heat  pipe  start-up  from 
frozen  state.  The  heat  pipe  under  consideration  was  rectangular  in 
geometry.  The  finite  element  method  was  used  to  solve  the  transient 
two-dimensional  solid  wall  and  frozen  liquid  regions.  A  one- dimensional 
analytical  approach  similar  to  that  of  Levy  (1968)  was  employed  to  solve 
the  steady- state  vapor  flow.  The  transient  temperature  variation  along  the 
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liquid- vapor  interface  was  presented  at  different  time  intervals.  The 
results  showed  that  it  would  take  about  150  seconds  for  a  heat  pipe  to  get 
to  the  normal  operation  state  from  the  frozen  state.  Since  the  1-D 
analytical  solution  for  vapor  flow  is  easy  to  use  and  gives  quite  accurate 
results,  it  was  a  good  approach  to  combine  this  method  together  with  the 
2- D  unsteady  solution  of  the  wall- liquid  regions.  This  is  because  the 
approaches  for  the  1-D  vapor  and  2-D  solid  regions  are  well  established, 
and  therefore  a  converged  solution  is  ensured. 

A  transient  heat  pipe  model  with  the  wall,  liquid- wick  and  vapor 
regions  was  proposed  by  Seo  and  El-Genk  (1988).  This  model  solves  the 
quasi- steady- state,  one- dimensional  mass,  momentum  and  energy  conservation 
equations  in  the  vapor  region  and  a  two-dimensional  transient  model  is  used 
in  the  wall  and  liquid- wick  regions.  This  model  assumes  a  compressible 
laminar  vapor  flow  and  an  incompressible  liquid  flow  in  the  porous  medium 
with  the  capillary  relationship  as  a  link  between  the  two  fluids,  but  the 
effect  of  porosity  in  the  liquid- wick  region  was  neglected.  Furthermore, 
the  energy  due  to  mass  injection  and  suction  to  and  from  the  vapor  region 
was  not  considered.  The  thermophysical  properties  were  assumed  to  be 
temperature  dependent  for  both  the  liquid  and  vapor.  The  solution  gives 
the  temperature,  pressure,  and  velocity  distributi  ns  in  the  heat  pipe,  as 
well  as  the  predictions  of  the  heat  pipe  operational  limits  during  both 
steady  and  transient  operations.  The  numerical  predictions  of  the 
temperatures  of  a  lithium  heat  pipe  with  a  radiative  boundary  condition  in 
the  condenser  were  in  good  agreement  with  experimental  data. 
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Issacci  et  al.  (1988)  studied  the  transient  behavior  of  the  vapor  flow 
in  a  rectangular  heat  pipe.  The  time  dependent  viscous  compressible 
elliptic  2-D  continuity,  momentum  and  energy  conservation  equations  were 
solved  numerically.  The  ideal  gas  law  is  used  to  model  the  compressibility 
of  the  vapor  and  the  pressure  and  temperature  are  linked  by  thermodynamic 
equilibrium  at  the  interface.  The  "SIMPLE"  numerical  scheme  developed  by 
Patankar  (1980)  was  used  for  the  computations.  Flow  reversal  was  detected 
in  the  condenser  and  adiabatic  regions  with  water  as  working  fluid,  but  no 
comparison  was  made  with  experimental  data. 

Faghri  (1988)  extended  the  solution  derived  by  Levy  (1968)  to  annular 
heat  pipes  by  using  a  similar  approach.  The  annular  heat  pipe  under 
consideration  has  two  concentric  pipes  of  unequal  diameters  attached  by  end 
caps  which  create  an  annular  vapor  space  between  the  two  pipes.  The  heat 
input  and  cooling  were  supplied  to  both  the  inner  and  outer  walls.  Similar 
results  as  obtained  by  Levy  (1968)  and  Bystrov  and  Popov  (1976)  were 
obtained  for  the  pressure  and  temperature  variation  along  the  pipe  at  the 
sonic  limit.  It  was  confirmed  that  the  friction  along  the  pipe  will 
accelerate  the  vapor  flow  and  make  the  sonic  limit  occur  at  a  smaller  heat 
input.  A  simple  equation  was  also  proposed  for  the  prediction  of  the  sonic 
limit  of  the  annular  heat  pipe.  In  addition,  the  two-dimensional 
numerical  parabolic  and  elliptic  solutions  were  also  obtained  for  the 
pressure  drop  along  the  pipe.  Based  on  the  observation  that  the  vapor  flow 
becomes  fully  developed  in  a  very  short  distance  from  the  evaporator  end 
cap,  three  analytical  similarity  expressions  were  obtained  for  the  pressure 
drop  in  the  three  different  sections  of  conventional  and  annular  heat 
pipes. 
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Experimental  Studies 

A  number  of  experiments  have  been  performed  both  on  actual  and 
simulated  heat  pipes  and  are  summarized  in  Table  1.2.  The  experimental 
data  on  the  actual  heat  pipes  are  usually  obtained  for  the  vapor  or  wall 
temperature  variations  along  the  heat  pipe.  The  pressure  profile  along  the 
pipe,  which  is  usually  difficult  to  obtain  for  actual  heat  pipes,  is 
measured  by  using  simulated  heat  pipes. 

Kemme  (1969),  Ivanovskii  et  al.  (1982)  and  Merrigan  et  al.  (1986) 
investigated  the  performance  of  liquid  metal  heat  pipes.  Kemme’s 
experiment  was  focused  on  the  situation  when  a  sodium  heat  pipe  is  working 
close  to  the  sonic  limit  (»  =  !)■  Ivanovskii  et  al.  studied  the  sodium 
heat  pipe  performance  at  different  working  temperatures  by  measuring  the 
vapor  saturation  temperature.  For  space  applications,  Merrigan  et  al. 
(1986)  performed  an  experiment  on  the  transient  performance  of  the  lithium 
heat  pipe  with  radiation  cooling  at  the  condenser. 

For  low  temperature  heat  pipes,  Ponnappan  and  Mahefkey  (1984),  Bianchi 
(1987),  Gernert  (1986)  and  Faghri  and  Thomas  (1988)  used  water  as  the 
working  fluid.  Ponnappan  and  Mahefkey  (1984)  studied  a  heat  pipe  with  a 
double-wall  artery  wick  structure  while  Bianchi  (1987)  studied  a 
conventional  heat  pipe  with  a  screen  mesh  wick.  Gernert’s  experiment  on 
the  heat  pipe  with  multiple  heat  sources  found  a  very  small  difference  in 
temperature  for  the  wall  and  vapor  temperatures  along  the  heat  pipe. 
Faghri  and  Thomas  (1988)  performed  an  experiment  on  an  annular  heat  pipe 
which  has  axial  grooves  on  the  inner  and  outer  pipe  walls  for  the  capillary 
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wick.  The  purpose  was  to  measure  the  capillary  limit  and  compare  the 
results  with  a  conventional  heat  pipe  with  the  same  outer  dimensions. 

The  most  widely  referenced  experiment  on  a  simulated  heat  pipe  is  that 
done  by  Quaile  and  Levy  (1972).  The  experiment  simulated  the  vapor  flow 
in  a  heat  pipe  condenser.  The  tube  was  closed  at  the  downstream  end  and 
the  fluid  was  removed  uniformly  by  suction  through  the  porous  wall.  They 
measured  the  axial  pressure  drop  in  a  porous  pipe  as  well  as  the  radial 
variations  of  the  axial  velocity  at  several  axial  locations.  The  pressure 
recovery  phenomenon  was  observed  in  that  experiment.  The  experimental  data 
from  the  pressure  profile  showed  a  very  good  agreement  with  the  theoretical 
solution  by  Veissberg  (1959)  and  Busse  (1967)  in  the  range  2.21  <  Rer  < 
5.0.  Recently,  Bowman  and  Hitchcock  (1988)  simulated  a  whole  heat  pipe  by 
employing  air  flow  with  injection  and  extraction  in  a  porous  pipe.  The 
results  were  given  for  the  pressure  and  velocity  variations  while  the  flow 
changed  from  laminar  to  turbulent. 
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1.3.  Mathematical  Formulation 

The  heat  pipe  model  under  consideration  (Fig.  1.1),  has  three  distinct 
regions  in  the  radial  direction  (wall,  liquid- wick  and  vapor)  as  well  as 
three  different  types  of  sections  in  the  axial  direction  (evaporator, 
adiabatic  and  condenser).  This  model  solves  the  two-dimensional 
conservation  of  mass,  momentum  and  energy  equations  for  the  heat  pipe  with 
single  and  multiple  heat  sources  under  the  following  assumptions: 

a.  The  compressible  vapor  flow  is  laminar  and  steady. 

b.  The  heat  transfer  through  the  liquid- wick  is  modeled  as  pure 
conduction  with  an  effective  thermal  conductivity. 

c.  The  properties  of  the  vapor,  the  liquid  and  the  solid  in  each 
region  are  constant  with  the  vapor  density  following  the  perfect 
gas  law. 

d.  Both  evaporation  and  condensation  are  considered  to  occur  at  the 
inner  radius  of  the  porous  medium. 

e.  At  the  vapor- liquid  interface,  the  vapor  is  at  its  thermodynamic 
equilibrium  temperature  corresponding  to  the  local  saturation 
vapor  .ressure. 

f.  At  the  liquid- vapor  and  wall- liquid  interfaces,  the  harmonic  mean 
of  the  thermal  conductivity  is  used  for  the  energy  equation. 

g.  The  flow  is  axisymmetric. 
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1.3.1.  Vapor  Region: 


The  conservation  of  mass  and  momentum  equations  given  by  Bird  et  al. 
(1960)  and  the  energy  equation  given  by  Bejian  (1984)  after  some 
simplification  are  as  follows: 


mass: 

$£(Vv)  +  ?  i^VV  =  0 


(i.i) 


r- component  momentum: 


(1.2) 


z- component  momentum: 


e/wv  _  u ry  n  d  x  d  "'y 
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energy  conservation: 
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where: 


The  viscous  dissipation,  <ji,  is  included  in  the  present  analysis.  V'e 
should  mention  here  that  in  equations  (1.2- 1.4),  the  terms  in  braces  {}  are 
associated  with  axial  diffusion  terms.  These  terms  are  neglected  when  the 
partially  parabolic  version  is  considered  but  are  included  in  the  elliptic 
version. 


Compressibility  of  Vapor: 

The  perfect  gas  law  is  employed  to  account  for  the  compressibility  of 
the  vapor.  The  vapor  density  is  obtained  from  the  following  equation: 

P 

Pv  =  (1.5a) 

v 

where  R  is  the  gas  constant. 


Since 


and 


we  have 
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7P 


v  =  c2 


where  c  is  the  velocity  of  sound  and  7  is  the  ratio  of  specific  heats. 
Usually  the  velocity  of  sound  changes  very  little  under  the  normal  working 
conditions  of  a  heat  pipe,  so  we  can  write  the  derivative  of  the  above 
equation  in  the  following  form: 


dp. 


p  d? 
A'v  y 


1 

r 


(1.5b) 


Equation  (1.5b)  is  used  in  the  density  correction  formula  p-p*+p' , 
where  p ’  is  the  density  correction  and  is  equal  to  ^P’,  P’  is  the  pressure 
correction,  p*  is  the  guessed  density  and  is  obtained  from  p  of  the 
previous  iteration.  As  long  as  a  converged  solution  can  be  obtained,  P’ 
approaches  zero,  therefore  the  approximate  p'  equation  (i.e.,  Eq.  (1.5b)) 
is  sufficient. 

1.3.2  Vail  and  Liquid- Vick  Regions 


In  these  regions,  only  heat  conduction  is  considered  and  the  energy 
conservation  equation  is  as  follows: 


1  d 
r  dr 


rk 


5T1 

3r 


8 

dz 


<7T' 

kdi 


=  0 


(1.6) 


where  the  thermal  conductivity  k  of  the  wall  is  different  from  that  of  the 
liquid- wick  structure. 

The  equation  for  the  effective  thermal  conductivity  proposed  by  Chi 
(1976)  is  used  here  for  the  liquid  saturated  screen  wicks,  and  is  as 
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follows: 


_  kitki  +  ks  '  U-OO^-kg)] 
keff  =  ki"+  ks  +  (l-e)(krkB) 


(1.7) 


In  the  above  equation,  the  e  is  the  wick  porosity  and  is  given  by  Chang 
(1987): 


e  =  1  - 


tAB 


where  A  =  d/w  and  B  =  d/t.  The  wire  diameter  is  d,  the  screen  thickness  is 
t,  and  the  opening  width  of  the  screen  is  w.  The  liquid  and  solid  wick 
thermal  conductivities  are  k-^  and  kg,  respectively. 


For  the  sintered  powder  wick,  the  following  equation  proposed  by  Dunn 
and  Reay  (1982)  is  used  to  calculate  the  effective  thermal  conductivity. 


'2  +  k1/kg  -  2c(l-k1/kg)' 
2  +  kx/ks  '+  e(l-k1/kg) 


(1.8) 


For  the  concentric  annulus  wick,  the  effective  thermal  conductivity  is 
obtained  from  the  following  equation  proposed  by  Dunn  and  Reay  (1982). 


=  k 


1 
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1.3.3.  The  Boundary  Conditions 


At  either  end  of  the  heat  pipe  (z  =  0,  z  =  L),  the  radial  and  axial 
velocities  and  the  axial  temperature  gradient  are  zero. 


w  =  v 


(1.9) 


The  pressure  is  fixed  to  zero  at  the  end  cap  of  the  evaporator  at  the 
liquid- vapor  interface  for  the  elliptic  solution  and  the  datum  pressure, 
PQ,  is  added  to  get  the  absolute  value. 


The  uniform  heat  flux  boundary  condition  is  specified  at  the  outer  pipe 
wall  surface  (r  =  R^)  of  the  evaporator  and  condenser  sections.  For  the 
evaporator  section  (0  <  z  <  L  )  or  each  of  the  evaporator  sections  in  the 

C 

case  of  multiple  heat  sources,  the  constant  heat  flux  boundary  condition  is 
as  follows: 

v 

'kw  3T  =  %,e 

For  the  adiabatic  section  (L  <  z  <  L  +  L  )  at  the  outer  pipe  wall 
surface  (r  =  Ry),  we  have 

^=°  (1.10) 


and  for  the  condenser  section  (L  +  L  <  z  <  L),  the  constant  heat  flux 

6  3- 

boundary  condition  is: 
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an  a 

,  w  e 

'kw  W~  =  qo,e  Tc 


At  the  wall  and  liquid- wick  interface  (r  =  Ry  +  <5^),  the  temperature 

and  heat  flux  at  both  sides  of  the  interface  should  be  the  same. 


Tt  =  T 
lw  w 


(1.11) 


91 


91 


=  k 


lw 


w  dr  eff  dr 


The  symmetry  boundary  conditions  are  applied  along  centerline  of  the 
pipe  (r=0): 


dv 


or 


\  =  w 


=  0 


(1.12) 


Along  the  liquid- vapor  interface  (r  =  Ry),  the  following  boundary 
conditions  are  specified:  no  slip  boundary  condition,  thermodynamic 
equilibrium,  energy  balance  heat  source,  and  blowing  and  suction 
velocities. 

For  the  no  slip  boundary  condition  at  the  liquid- vapor  interface 


w  =  0 
v 


(1.13) 


For  thermodynamic  equilibrium,  the  Clapeyron  equation  is  used  in  the 
following  form: 


dP  dT  h, 
v  v  f 
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or  it  can  be  integrated  in  terms  of  Tq  and  Pq,  which  results  in 

1 

T 


T  = 
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(1.14) 
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The  above  equation  is  used  to  determine  the  vapor  temperature  at  the 
vapor- liquid  interface. 


Making  an  energy  balance  at  the  interface  results  in 
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^  >  0  Suction  (condensatiun) 

ri  “  jcT  {  v^  <  0  Blowing  (evaporation) 


(1.15) 


(1.16) 


(1.17) 


A  summary  of  the  boundary  conditions  are  also  listed  in  Table  1.3. 
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Table  1.3  Boundary  Conditions  for  Heat  Pipe  Analysis 


evaporator 
(0  <  z  <  Le) 

adiabatic 
(Le  <  2  <  te  +  La) 

condenser 
(Le  +  La  <  z  <  L) 

outer  pipe 
vail 

(r  =  »„) 

a 

f 

a> 

O 

11 

w  u 

*e 

"kv  Tr  =  %,e  ^ 

wall  and 
liquid- wick 
interface 

(r  =  »v  -  «!> 

«lw 

Tl«  =  V  kvT?  1  keff  TT 

liquid- vapor 
interface 

(r  =  »,) 

t  «1. 

“y  =  °>  'i  =  j-y’  1j  =  W,fg  =  kv~»  -  kefOT  > 

1 

Ti”  <  D  P_ 

1_  .  R  In  Z 

T0  ®fg  P0 

centerline  of 
the  pipe 
(r  -  0) 

dv 

v  &l 

=  '  -5?  =  ° 

both  ends  of 

f 

the  pipe 
(z  =  0,  L) 

i 

• *  v  =  i =  0 
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1.3.4.  The  Dimensionless  Variables 


The  solution  for  the  parametric  study  will  be  presented  for  the 
following  dimensionless  dependent  variables: 


+ 

w 


2(P  -  P0) 


where  p ^  is  obtained  at  Tq  and  Pq.  The  mean  heat  flux  at  the  liquid- vapor 
interface  without  considering  the  axial  heat  conduction  in  the  wall  and 
liquid-wick  regiouo  is  qm,  which  is  specified  as  follows: 

at  the  evaporator  and  adiabatic  sections, 

qm,e  "  %,e 

and  at  the  condenser  section, 

qm,c  ”  %,c 

where  q  and  q  are  the  heat  fluxes  specified  at  the  outer  pipe  wall  of 

U  J  V  V  J  V 

the  evaporator  and  condenser  sections.  The  heat  flux  across  the 
liquid- vapor  interface,  q^,  is  calculated  from  equation  (1.15). 

These  results  are  presented  in  terms  of  the  dimensionless  independent 
variables  and  parameters  as  follows: 
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The  dimensionless  parameters  proposed  here  are  based  on  the  conjugate 
heat  transfer  analysis  of  flow  in  a  pipe  (Section  II)  and  the  annular  heat 
pipe  analysis  by  Faghri  and  Parvani  (1988). 
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1.4  Solution  Procedure 


The  heat  pipe  problem  is  solved  as  a  convection- conduction  problem 
throughout  the  entire  domain  by  solving  one  generalized  energy  equation 
with  different  thermal  diffusion  coefficients.  The  velocity  in  the  solid 
wall  and  the  liquid- wick  matrix  is  specified  to  be  zero  so  that  the 
analysis  in  these  two  regions  becomes  purely  a  conduction  problem.  Since 
the  liquid  velocity  in  the  porous  medium  is  much  slower  compared  to  the 
vapor  flow,  the  zero  vapor  velocity  boundary  condition  at  the  liquid- vapor 
interface  and  neglecting  the  convective  term  in  the  energy  equation  of  the 
liquid- wick  region  should  not  cause  a  large  accuracy  problem.  The  validity 
of  the  results  will  be  checked  with  experimental  data. 

The  generalized  PHOENICS  computer  code  is  used  in  the  present  analysis, 
which  employs  the  finite-difference  iterative  method  of  solution  developed 
by  Spalding  et  al  (1980).  The  elliptic  solutions  of  the  mass,  momentum 
and  energy  conservation  equations  with  boundary  conditions  given  in 
equations  (1.9-1.17)  were  obtained.  The  partially  parabolic  solution  was 
also  obtained  by  neglecting  the  axial  diffusion  terms  in  the  momentum  and 
energy  equations.  The  solution  procedure  is  based  on  a  line-by-line 
iteration  method  in  the  axial  direction  and  the  Jacobi  point- by- point 
procedure  in  the  radial  direction.  The  "SIMPLEST"  method  is  employed  for 
the  momentum  equations,  in  which  the  finite- domain  coefficients  contain 
only  diffusion  contributions,  and  where  the  convection  terms  are  added  to 
the  linearised  source  term  of  the  equations. 

The  pressure  field  is  solved  by  the  whole- field  pressure  correction 
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algorithm  derived  by  Markatos  et  al.  (1982).  The  pressure  correction 
equation  is  deduced  from  the  finite- domain  form  of  the  continuity  equation, 
which  assumes  that  the  flow  coming  into  the  system  is  due  to  the  pressure 
difference.  The  pressure  field  is  first  assumed  and  then  the  velocities 
and  temperatures  are  solved.  This  step  is  repeated  line-by-line  to  the 
completion  of  a  sweep.  The  pressure  correction  equation  is  then  solved 
using  the  mass  errors  that  have  been  calculated  during  the  sweep  and  other 
variables  are  updated  accordingly.  A  new  sweep  will  start  until 
convergence  is  attained.  In  this  numerical  scheme,  once  the  velocity 
boundary  condition  is  specified  no  pressure  boundary  condition  is  needed. 
For  compressible  vapor  flow,  the  mass  and  momentum  equations  are  coupled 
with  the  energy  equation  by  the  perfect  gas  law.  These  four  equations  are 
solved  sequentially  by  iteration. 


Since  there  is  a  change  of  phase  at  the  liquid- vapor  interface,  the 
energy  equation  is  no  longer  continuous  due  to  the  latent  heat  of 
evaporation  or  condensation.  To  make  an  energy  balance,  we  can  include  the 


term  mh^  as  a  heat  sink  at  the  liquid- vapor  interface  in  the  evaporator 
and  as  a  heat  source  at  the  interface  in  the  condenser  section.  Therefore, 


the  sign  of  qi  in  the  evaporator  should  be  negative  and  positive  in  the 
condenser. 


The  governing  equation  of  the  vapor  flow  is  first  solved  by  assuming 
that  the  heat  flux  is  uniform  at  the  liquid- vapor  interface  based  on  the 
total  heat  input  at  the  outer  wall.  The  rate  of  evaporation,  condensation 


and  the  velocity  are  then  calculated  by  m  =  q^/h^  and  v^ 


These 


values  are  used  in  the  vapor  mass  and  momentum  conservation  equations. 
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Once  the  converging  solution  is  obtained  after  a  few  iterations,  a  new  q- 
is  calculated  through  equation  (1.15)  and  used  for  further  iterations.  At 
the  same  time,  the  thermodynamic  equilibrium  is  checked  and  corrected 
during  each  iteration.  The  iteration  will  continue  until  fully  converged 
results  are  obtained.  The  pressure  at  the  liquid- vapor  interface  at  the 
end  of  the  evaporator,  Pq,  is  taken  as  the  datura  pressure  and  does  not 
change.  The  corresponding  saturation  temperature  Tq  is  assumed  to  be  the 
initial  temperature  for  the  whole  calculation  domain. 

The  accuracy  of  the  numerical  solution  is  checked  with  experimental 
data  and  the  convergence  is  assured  in  two  ways: 

1)  The  sum  of  the  absolute  value  of  the  residuals  should  decrease  as  the 
sweep  number  increases. 

2)  The  spot  value  should  approach  a  constant  value  as  the  sweep  number 
increases. 

For  each  grid  point,  a  residual  is  defined  as  the  error  that  occurs 
when  the  current  values  of  the  dependent  variables  are  inserted  into  the 
discretization  equation.  The  error  is  due  to  the  fact  that  the  current 
values  of  the  dependent  variables  do  not  satisfy  the  discretization 
equation  exactly.  The  number  of  sweeps  changes  from  150  to  350  based  on 
each  different  case.  The  spot  value  was  used  to  monitor  the  change  of  each 
variable  with  the  increase  of  the  sweep  number  at  one  particular  location 
in  the  domain.  The  location  of  the  spot  value  was  chosen  in  the  vapor 
region  of  the  condenser  section  next  to  the  adiabatic  section. 


As  a  common  approach,  a  course  grid  size  is  first  chosen  to  test  the 
program,  and  a  fine  grid  spacing  is  employed  for  the  final  solution.  A 
uniform  grid  size  is  used  for  the  axial  direction,  and  three  different 
uniform  grid  sizes  are  used  for  the  vapor,  liquid- wick,  and  wall  regions. 
For  the  numerical  results  presented  for  the  parametric  study,  the  grid 
sizes  for  the  cases  presented  are  chosen  as  follows: 


30  x  20  x  10  x  10  =  (axial)  x  (radial  vapor)  x  (radial  liquid- wick)  x 
(radial  solid  wall). 


1.5  Model  Verification  Versus  Experimental  Data 


To  verify  the  numerical  predictions,  the  results  are  compared  to  four 
cases  of  existing  experimental  data  reported  by  several  investigators  as 
listed  in  Table  1.4.  The  related  experimental  specifications  and 
properties  for  the  numerical  computations  are  also  listed  in  Tables  1.4  and 
1.5,  respectively. 

The  numerical  model  was  first  compared  to  the  experimental  data 

reported  by  Ivanovskii  et  al.  (1982)  for  a  cylindrical  sodium  heat  pipe 

(Case  1).  The  heat  pipe  was  provided  with  a  compound  wick  of  the  type  with 

a  ring-shaped  gap  for  the  flow  of  liquid.  The  effective  thermal 

conductivity  for  the  computation  is  chosen  to  be  equal  to  the  liquid 

thermal  conductivity  proposed  by  Dunn  and  Reay  (1982).  The  method  of 

measuring  the  temperature  distribution  was  to  place  a  movable 

micro- thermocouple  directly  in  the  vapor  channel.  The  thermocouple  was 

provided  with  a  special  capillary  device  to  keep  it  wetted  by  the 

condensate,  so  the  temperature  readings  correspond  to  the  saturation 

temperature  which  matches  the  numerical  model  at  the  interface.  For  the 

numerical  computation,  the  constant  heat  flux  with  a  total  heat  input  of 

560  V  is  specified  at  the  outer  wall  of  the  evaporator.  The  pipe  wall  is 

assumed  to  be  made  of  stainless  steel  and  has  a  thickness  of  6  =1  mm.  At 

w 

the  condenser,  a  constant  heat  flux  based  on  the  total  heat  input  at  the 
outer  wall  of  the  evaporator  section  is  also  specified.  Fig.  1.2  compares 
the  numerical  results  and  the  experimental  data  for  the  vapor  saturation 
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TABLE  1.4  Experimental  Heat  Pipe  Specifications 


Case  No. 

1 

2 

3 

4 

References 

Ivanovskii 

Ivanovskii 

Kemme 

Gernert 

et  al.  (1982) 

et  al.  (1982) 

(1969) 

(1986) 

Working 

Fluid 

Sodium 

Sodium 

Sodium 

Water 

M») 

0.1 

0.1 

0.143 

0.1,  0.1 

M"* 

0.05 

0.05 

0.06 

0.4,  0.3 

MB> 

0.35 

0.55 

1.08 

0.6 

Vapor  Channel 
Radius  (mm) 

7.0 

7.0 

5.7 

11.0 

Vick  Type  Gap( 

Ring- Shaped 
screen  mesh) 

Ring- Shaped 
Gap(screen  mesh) 

Screen  Mesh 

Sintered 

Powder 

Mesh  Number 
(per  in2) 

not  reported 

not  reported 

400 

325 

Vick  Thickness 
(tfpinm) 

0.5 

0.5 

0.15 

0.76 

Wall  Stainless 

Material  Steel(assumed) 

Stainless 

Steel(assumed) 

Stainless 

Steel 

Copper 

Wall  Thickness 
(<5w,mm) 

1.0 

(assumed) 

1.0 

(assumed) 

0.9 

1.6 

Total  Heat 
Input  (Q , V) 

560 

1000 

6400 

200  +  800 
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TABLE  1.5  Experimental  Heat  Pipe  Properties 


Case  No. 

1 

2 

3 

4 

Properties 

y°c) 

545 

583 

692 

94 

p0C/"!) 

1300 

2476 

12,160 

85,710 

hf6xlO_6(J/kg) 

4.22 

4.182 

4.08 

2.27 

/iTxl02kg/n3) 

0.461 

0.843 

3.857 

51.06 

/ivxl02N  s/m2) 

194.9 

191.43 

192.89 

121.42 

c„  (J/kg'K) 

Fv 

2583.0 

2648.0 

2710.0 

2010.1 

t»S 

i 

s 

:* 

to 

19.0 

19.0 

19.0 

394.0 

k1(W/B  K) 

66.18 

66.18 

59.54 

0.682 

ke££(W/m-K) 

66.18 

66.18 

45.45 

240.2 

kv(V/m-K) 

0.0352 

0.0353 

0.044 

0.0244 

e 

0.33 

0.33 

0.74 

0.30 

grid  number 

(radial  x  axial) 

35  x  50 

35  x  70 

35  x  50 

40  x  60 
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temperature  along  the  heat  pipe  for  Case  1.  These  particular  experimental 
data  were  chosen  for  comparison  because  they  were  more  readable  than  other 
data  presented  by  Ivanovskii  et  al.  (1982).  As  the  results  show,  the 
present  compressible  elliptic  and  partially  parabolic  models  give  accurate 
predictions  of  the  temperature  profile  compared  to  the  experimental  data 
with  a  maximum  deviation  of  3°C.  According  to  Ivanovskii  et  al.  (1982), 
the  experimental  data  for  the  heat  transfer  rate  were  measured  with  an 
accuracy  of  6  to  107.,  so  the  deviations  of  the  present  compressible  models 
are  within  the  range  of  experimental  accuracy.  For  the  incompressible 
model,  however,  there  is  a  maximum  deviation  of  about  6°C  at  the  inlet  of 
the  condenser  section.  Therefore,  we  need  to  include  the  effect  of  the 
compressibility  of  the  vapor  in  the  analysis. 

Figure  1.3  shows  the  variation  of  the  pressure  drop  distribution  along 
the  liquid- vapor  interface  for  Case  1.  The  pressure  drop  reaches  its 
maximum  value  at  the  exit  of  the  adiabatic  section  and  then  recovers  about 
557.  in  the  condenser.  Since  the  Clapeyron  equation  is  the  link  between  the 
temperature  and  the  pressure  at  the  liquid- vapor  interface,  the  pressure 
profile  is  similar  to  the  temperature  distribution.  The  trend  of  the 
pressure  profile  also  agrees  with  the  existing  numerical  results  obtained 
by  Tien  and  Rohani  (1974). 
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P-P0,  N/m 
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Fig.  1.4  shows  the  variation  of  the  Mach  number  along  the  centerline  of 
the  pipe  with  a  maximum  value  of  0.6  at  the  inlet  of  the  condenser  for  the 
compressible  elliptic  model.  The  trend  is  generally  in  agreement  with  the 
results  obtained  by  Levy  (1968)  for  a  cylindrical  heat  pipe  and  Faghri 
(1988)  for  an  annular  heat  pipe  with  a  1-D  compressible  model. 

Figures  (1.5- 1.7)  shows  the  numerical  results  for  the  axial 
temperature,  pressure  and  Mach  number  variations,  respectively,  which 
correspond  to  Case  2  of  the  experimental  data  by  Ivanovskii  (1982).  For 
this  case,  the  heat  input  was  Q=1000  V  and  other  geometric  parameters  are 
indicated  in  Table  1.4.  Since  the  thickness  of  the  wall  and  the  type  of 
materials  were  not  mentioned  in  the  reference,  numerical  calculations  were 
made  with  two  different  thermal  conductivities  and  thicknesses  for  the  pipe 
wall  which  are  believed  to  be  commonly  used  by  heat  pipe  manufacturers. 
The  results  showed  that  the  effect  of  axial  conduction  inside  the  wall  and 
liquid- wick  is  not  important  in  this  range,  so  only  the  results  for  the 
thermal  conductivity  and  wall  thickness  reported  in  Tables  1.4  and  1.5  are 
given.  The  results  of  the  compressible  and  incompressible  models  for  the 
vapor  flow  are  presented  versus  the  experimental  data  for  the  interface 
temperature  in  Fig.  1.5.  It  shows  that  both  models  give  a  very  good 
prediction  in  the  evaporator  region,  but  the  incompressible  model 
overpredicts  the  data  in  the  adiabatic  and  condenser  regions  while  the 
compressible  model  underpredicts  the  data  in  the  start  of  the  condenser 
region  and  overpredicts  it  near  the  end.  In  general,  the  compressible 
model  gives  a  better  prediction.  Fig.  1.6  shows  the  pressure  variation 
along  the  heat  pipe  for  the  compressible  and  incompressible  models  for  Case 
2.  The  incompressible  model  gives  a  maximum  deviation  of  about  187. 
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compressible,  elliptic  compr^ssjble^parabolic 


number  along  the  centerline  of  the  sodium  heat  pipe 


■  exp.  data  by  Ivanovskii  (1982)  incompressible,  dliptic 
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Fig.  1.5  The  axial  interface  temperature  profile  along  the  sodium  heat  pipe 
with  Q=1000  W,  Rv=0.007  m,  Le=0.1  m,  L,=0.05  m,Lc=0.55  m 
kj=66.2  W/m2-K,  ks=19.0  W/m  -K,  <5j=0.0005  m,  <5W=0.001  m 


compressible,  e 
incompressible, 


number  profile  along  centerline  of  the  sodium  heat  pipe 
W.  =0.007  m.  L«=0.1  m,  L^=0.05  m,  L~=0.55  m 


compared  to  the  compressible  model.  Fig.  1.7  presents  the  Mach  number 
variation  obtained  from  the  compressible  model  with  the  maximum  Mach  number 
of  0.54. 

Figures  1.8  and  1.9  show  the  numerical  results  for  the  axial 
temperature  and  Mach  number  variations,  respectively,  which  correspond  to 
Case  3  of  the  experimental  data  by  Kemme  (1969).  The  heat  pipe  was  about 
1.3  m  long  and  5.7  mm  I.D.  with  a  screen  wrap  type  wick  of  thickness  = 
0.15  mm  and  stainless  steel  wall  of  thickness  6 =  0.9  mm.  Heat  was  added 
to  the  evaporator  section  of  the  heat  pipe  with  an  induction  coil  while  the 
heat  was  removed  from  the  condenser  section  by  conduction  through  a  gas  gap 
to  a  water  calorimeter.  The  effective  thermal  conductivity  of  the  wick  is 
calculated  from  Eq.  (1.7).  For  this  case,  three  different  sets  of 
experimental  data  of  the  outer  wall  temperature  were  obtained  for  subsonic, 
sonic  and  supersonic  vapor  flow  in  the  sodium  heat  pipe.  During  the 
experiment,  the  heat  input  was  fixed  to  6.4  kV  and  the  working  temperature 
was  decreased  by  changing  the  cooling  condition  so  the  choking  condition 
could  be  reached. 

For  the  present  numerical  analysis,  the  temperature  at  the  end  cap  of 
the  evaporator  is  fixed  to  the  experimental  value  of  692°C,  and  one 
steady- state  solution  is  obtained  with  a  Mach  number  of  1.0  at  the  exit  of 
the  adiabatic  section  as  shown  in  Fig.  1.9  for  the  compressible  model.  The 
maximum  Mach  number  of  1.09  was  found  in  the  inlet  region  of  the  condenser. 
Fig.  1.8  presents  the  experimental  data  of  the  wall  temperature  and  the 
numerical  solutions  of  the  compressible  and  incompressible  models.  The 
present  results  of  both  models  are  in  general  agreement  with  Kemme’ s  data 
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compressible,  elliptic 


number  profile  along  the  sodium  heat  pipe 
0057  m,  U=0.143  m,  L.^0.06  m,  Lc=1.08  m 
=19.0  W/m -K  <5i=O.OOOl5  m,  0.0009  m 


for  the  sonic  limit  case  except  that  the  condenser  wall  temperature  is 
higher  than  the  experimental  data.  This  is  probably  because  the  condenser 
cooling  may  affect  the  accuracy  of  the  thermocouple  reading  and  makes  the 
reading  lower  than  the  actual  value.  However,  according  to  the  comparison 
it  seems  that  the  present  model  can  still  predict  the  general  trend  of  the 
heat  pipe  performance  even  if  the  Mach  number  of  the  vapor  flow  is  high. 

To  further  check  the  validity  of  the  present  numerical  analysis,  the 
code  was  also  modified  to  predict  the  performance  of  the  water  heat  pipe 
with  multiple  heat  sources.  The  results  were  compared  with  the 
experimental  data  obtained  by  Gernert  (1986).  The  experimental  heat  pipe 
under  consideration  had  two  evaporators  and  one  condenser  as  illustrated  in 
Fig.  1.10.  The  heat  was  provided  at  the  evaporator  by  two  electric  heater 
blocks.  The  condenser  section  was  fitted  with  a  water  cooled  jacket.  The 
heat  pipe  had  a  sintered  powder  wick  and  the  effective  thermal  conductivity 
was  calculated  from  the  equation  [Eq.(1.8)]  proposed  by  Dunn  and  Reay 
(1982).  A  multipoint  thermocouple  was  installed  along  the  centerline  of 
the  vapor  space  to  measure  the  variations  of  the  vapor  temperature.  To 
measure  the  outer  wall  temperature,  the  outer  wall  of  the  evaporator  and 
condenser  sections  were  grooved,  and  thermocouples  were  soldered  in  the 
grooves  to  the  pipe  wall.  According  to  the  experimental  conditions,  two 
constant  heat  fluxes  with  heat  inputs  of  200  V  and  800  V  were  specified  at 
the  two  evaporator  sections.  For  the  condenser  section,  a  constant  heat 
flux  based  on  the  total  heat  input  is  also  specified  as  a  boundary 
condition  for  the  numerical  computation. 
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Fig.  1.10  The  axial  temperature  profile  along  the  water  heat  pipe 

with  multiple  heat  sources 


Figure  1.10  shows  the  numerical  results  and  the  experimental  data  for 
the  outer  wall  and  the  vapor  temperature  of  the  water  heat  pipe  with 
multiple  heat  sources  for  Case  4.  The  present  model  gives  a  very  uniform 
vapor  temperature  and  a  very  small  outer  wall  temperature  variation.  The 
measured  vapor  temperature  is  also  fairly  uniform  with  less  than  a  1°C 
temperature  difference  along  the  pipe.  The  use  of  copper  as  the  wall  and 
wick  materials  and  the  sintered  powder  wick  give  a  very  good  thermal 
conductivity,  thus  the  radial  temperature  drop  is  very  small  through  the 
wall  and  liquid- wick. 

For  the  outer  wall  temperature,  two  experimental  data  points  found  in 
the  adiabatic  section  and  one  in  the  800  V  evaporator  are  fairly  close  to 
the  numerical  solutions,  but  the  temperature  readings  in  the  condenser 
section  are  almost  6°C  below  the  present  predictions.  Since  the  length  of 
the  condenser  is  much  longer  than  that  of  the  total  length  of  the 
evaporator,  the  heat  flux  through  the  pipe  wall  of  the  condenser  is  smaller 
than  that  in  the  evaporator,  which  means  that  the  heat  flux  through  the 
pipe  wall  of  the  condenser  is  smaller  than  that  in  the  evaporator  with  800 
V.  Thus,  it  is  expected  to  have  a  smaller  temperature  gradient  through  the 
pipe  wall  in  the  condenser  section.  It  seems  the  measurements  were  not 
very  accurate  in  the  condenser  section  probably  due  to  the  influence  of  the 
cooling  water  on  the  thermocouple  readings  which  resulted  in  the  readings 
being  lower  than  what  is  expected.  Usually  it  is  very  difficult  to  measure 
the  condenser  and  evaporator  wall  temperature  because  the  cooling  water  may 
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lower  the  actual  measurement  while  the  heating  may  cause  higner  temperature 
readings.  Furthermore,  the  accuracy  of  the  thermocouple  reading  may  also 
cause  problems  when  the  readings  differ  by  only  1  or  2°C. 

Figure  1.11  shows  the  pressure  variation  along  the  vapor- liquid 
interface  of  the  water  heat  pipe  for  Case  4.  A  very  small  pressure  drop  is 
noticed  along  the  200  V  evaporator  and  the  adjacent  adiabatic  section,  but 
in  the  800  V  evaporator,  the  additional  amount  of  heat  flux  causes  more 
liquid  to  be  evaporated  into  the  vapor  channel  thus  leading  to  a 
significant  pressure  drop  along  this  section.  In  the  condenser  section, 
about  807,  of  the  pressure  drop  is  recovered.  As  the  results  show,  the 
total  pressure  drop  is  very  small  which  means  that  a  very  uniform  vapor 
temperature  profile  can  be  obtained.  The  uniform  temperature  profile  is  a 
result  of  the  thermodynamic  equilibrium  between  the  pressure  and 
temperature  at  the  liquid- vapor  interface  coupled  with  the  fact  that  water 
has  a  high  static  vapor  pressure  and  vapor  density  under  normal  working 
conditions.  The  phenomenon  observed  here  is  in  general  agreement  with  what 
is  observed  for  a  heat  pipe  with  one  evaporator  section.  Ve  believe  that 
the  present  model  can  predict  the  general  performance  of  heat  pipes  with 
multiple  heat  sources  and  will  also  provide  a  guideline  for  future 
experiments  in  this  respect. 
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1.6  Parametric  Analysis 


The  heat  pipe  under  consideration  for  analysis  in  this  section  has  an 
evaporator,  adiabatic  and  condenser  section  each  with  a  corresponding 
length  of  :  L  =  0.2  m,  L  =  0.1  m,  L  =  0.3  m.  The  dimensions  of  the  pipe 
are  vapor  core  radius  Ry  =  8.6  mm,  wall  thickness  6  =  2  mm,  liquid  wick 
thickness  6-^  =  1.44  mm,  total  length  L  =  0.6  m.  The  wall  and  wick 
materials  are  assumed  to  be  stainless  steel  with  the  screen  wrap  wick  of 
porosity  e  =  0.7.  The  cases  presented  in  this  section  are  listed  in  Table 
1.6,  and  the  corresponding  properties  are  listed  in  Table  1.7.  The 
properties  are  obtained  at  temperature  Tq  and  are  assumed  to  be  constant 
along  the  pipe.  Since  the  radial  Reynolds  number  along  vapor- liquid 
interface  changes,  the  average  radial  Reynolds  numbers  are  calculated  based 
on  the  total  heat  input  at  the  outer  wall  of  the  heat  pipe  for  the  purpose 
of  comparison,  which  are  reported  in  Table  1.6.  The  heat  input  in  the 
evaporator  is  chosen  so  that  the  Mach  number  has  a  wide  range  for  the 
sodium  heat  pipe,  and  the  total  heat  input  is  smaller  than  the  heat  pipe 
limits. 
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Table  1.6  Case  Specifications  for  Parametric  Analysis 


Case  no. 

Working 

Fluid 

Re„ 

a 

Re  c 
r,e 

Rer  , 
r,c 

m 

T0(°C) 

P0(N/„2) 

A 

sodium 

186 

-2 

1.33 

213.4 

535 

1023 

B 

sodium 

372.1 

-4 

2.67 

426.7 

535 

1023 

C 

sodium 

558.1 

-6 

4.0 

640.2 

535 

1023 

D 

water 

4419.0 

-50 

33.3 

1647 

100 

101300 

Table  1.7  Properties  for  Parametric  Analysis 


1.  Sodium  Heat  Pipe  (at  Tn  =  535°C) 


hr  xlO  ^ 

fg 

(J/kg) 

4.237 

/?vxl02 

(kg/m3) 

3.713 

v10? 

(N-s/m2) 

207.5 

(J/lg-K) 

2555 

ks 

(W/m-K) 

17.4 

kl 

(W/m-K) 

59.5 

k  xlO2 

V 

(W/m-K) 

3.53 

2. 

Water  Heat  Pipe  (at  Tq 

i  =  100°C) 

hr  xlO  ® 

fg 

pvxl02 

fi  xlO^ 

'v 

CPv 

ks 

kl 

k  xlO2 

V 

(J/kg) 

(kg/m3) 

(N-s/m2) 

(J/kg-K) 

(W/m-K) 

(W/m-K) 

(W/m-K) 

2.257 

59.77 

120.3 

2028 

19 

.685 

2.479 

1.6.1.  Effect  of  Conjugate  Axial  Conduction 

Table  1.8  shows  the  six  cases  which  have  been  investigated  for  the 
conjugate  heat  transfer  effect  using  the  vapor  properties  of  sodium.  The 
property  and  parameter  specifications  needed  for  computation  are  the  same 
as  that  of  Case  B  (Table  1.6)  except  for  the  thermal  conductivities  and  the 
wall  and  wick  thicknesses.  For  all  of  the  cases,  the  vapor  thermal 
conductivity  and  vapor  core  radius  are  held  constant,  which  are  =  0.0353 
(V/m-K)  and  Ry  =  0.0086  m. 


Table  1.8  Case  Specifications  for  Conjugate  Effect  Analysis 


Case  no. 

“iv 

Kwl 

A 

w 

Case  No. 

Klv 

*»1 

41 

A 

w 

B 

983. 

0.69 

0.084 

0.116 

B3 

100 

1000 

0.1 

0.1 

Bl 

100 

10 

0.1 

0.1 

B4 

100 

10 

0.2 

0.2 

B2 

100 

100 

0.1 

0.1 

B5 

100 

10 

0.01 

0.01 

Figures  1.12  and  1.13  show  the  effect  of  the  thermal  conductivity  ratio 
on  the  performance  of  a  sodium  heat  pipe  with  the  pipe  wall  and 

liquid-wick  dimensions  held  at  A^  =  A^  =  0.1,  Re^  g  =  -4.0,  and  =  100. 

Fig.  1.12  presents  the  interfacial  heat  flux  variation  along  the 

liquid- vapor  interface  with  three  different  thermal  conductivity  ratios: 

=  10,  100,  1000  (Case  Bl,  B2,  B3).  In  addition,  the  solution  for  Case 
B  is  also  presented.  Fig.  1.13  shows  the  outer  wall  and  interface 
temperature  variations  for  the  same  cases.  The  effect  of  conjugate  axial 
conduction  becomes  very  significant  as  increases.  When  K  ^  =  10,  the 

results  show  a  fairly  uniform  distribution  of  heat  flux  with  the  ratio 
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}f  the  thermal  conductivity  ratio  Kwi  on  the  interfacial 
variations  along  the  heat  pipe  with  Rere=-4.0 


thermal  conductivity  ratio  K  on  the  axial  temperature 


close  to  unity  at  both  the  evaporator  and  condenser,  and  zero  in  the 
adiabatic  section.  As  increases,  the  results  show  that  the  conjugate 
effect  makes  the  heat  flux  less  uniform  and  the  heat  flux  deviates  from 
unity.  The  axial  wall  conduction  in  the  evaporator  unids  to  transfer  the 
heat  to  the  adiabatic  section.  Since  the  constant  heat  flux  boundary 
condition  is  specified  both  at  the  outer  wall  of  the  evaporator  and 
condenser,  the  interfacial  heat  flux  in  the  evaporator  becomes  less  than  1. 
In  the  condenser  section,  the  heat  is  conducted  from  the  evaporator,  so 
that  less  heat  flux  is  needed  at  the  interface  to  satisfy  the  outer  wall 
constant  heat  flux  boundary  condition.  As  the  results  show,  the 
interfacial  heat  flux  starts  from  a  negative  number  at  the  end  of  the 
adiabatic  section  and  gradually  increases  to  unity  at  the  end  of  the 
condenser.  At  =  1000,  a  significant  variation  of  the  interfacial  heat 
flux  can  be  noticed  due  to  axial  conduction.  Obviously,  the  inclusion  of 
axial  wall  and  liquid  wick  conduction  is  very  important  for  this  case. 
Also,  the  results  for  Case  B  are  almost  identical  to  this  solution  when  no 
axial  conduction  is  considered.  Thus,  the  conjugate  effect  of  this  case 
can  be  excluded  while  discussing  other  effects. 

The  interfacial  heat  flux  should  be  zero  at  the  adiabatic  section  if 
there  is  no  axial  conduction,  but  because  of  the  inclusion  of  axial 
conduction,  the  heat  flux  along  the  adiabatic  section  is  no  longer  zero. 
Since  the  local  qm  at  the  adiabatic  section  is  zero,  the  q  of  the 
evaporator  is  used  here  for  the  calculation  of  q - /q  in  the  adiabatic 
section.  The  negative  value  of  q^/q  in  some  regions  of  the  adiabatic 
section  is  because  the  sign  of  q.  changes  in  the  adiabatic  section. 
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Actually,  because  of  the  axial  conduction  in  the  pipe,  the  region  of  the 
adiabatic  section  near  the  condenser  acts  as  part  of  the  condenser, 
likewise  the  region  adjacent  to  the  evaporator  acts  as  part  of  the 
evaporator.  Only  one  point  in  the  adiabatic  section  is  truly  adiabatic 

For  =  1000,  the  interface  temperature  distribution  becomes  more 
uniform  and  correspondingly,  the  outer  wall  temperature  varies  more 
smoothly  as  shown  in  Fig.  1.13.  Therefore,  it  is  always  beneficial  to  have 
a  large  for  better  heat  pipe  performance. 

Figures  1.14  and  1.15  show  the  effect  of  the  wall  and  liquid- wick 

thicknesses  on  the  heat  pipe  performance.  Three  different  thicknesses  A^  = 

A^  =  0.2,  0.1,  0.01  (Case  B4,  Bl,  B5)  were  chosen  with  =  100,  =  10 

and  Re  =  -4.0.  Fig.  1.14  presents  the  heat  flux  variation  along  the 
r ,  6 

liquid- vapor  interface.  When  the  thickness  is  very  small  (Ay  =  A^  =  0.01), 
the  results  show  that  the  heat  flux  at  both  the  evaporator  and  condenser 
almost  equals  unity  which  means  the  axial  conduction  is  insignificant.  As 
A  increases,  the  effect  of  axial  conduction  becomes  more  pronounced.  Fig. 
1.15  gives  the  interface  and  outer  wall  temperature  distribution.  The 
outer  wall  temperature  should  be  equal  to  the  interface  temperature  at  the 
adiabatic  section  if  there  is  no  axi'-  .uction,  but  as  the  results  show, 
there  is  a  temperature  difference  for  the  cases  presented  and  the  case  Ay  = 
Aj  =  0.2  has  the  largest  difference.  However,  no  significant  effect  of 
axial  conduction  is  noticed  within  the  range  presented.  It  is  confirmed 
here  that  for  small  values  of  A,  the  axial  heat  conduction  through  the  wall 
and  liquid- wick  can  be  neglected  because  the  results  approach  that  of  a 
constant  heat  flux  at  the  interface.  It  is  also  interesting  to  note  that 
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thickness  on  the  interfacial  heat  flux 


thickness  on  the  axial  temperature 


the  outer  wall  temperature  in  the  middle  of  the  adiabatic  section  is  fairly 
close  to  the  vapor  temperature,  which  is  a  common  practice  to  estimate  the 
vapor  temperature  by  measuring  the  outer  wall  temperature  in  the  middle  of 
the  adiabatic  section. 

1.6.2  Effect  of  Vapor  Compressibility 

Figures  1.16  through  1.18  show  the  effect  of  compressibility  on  the 

sodium  heat  pipe  vapor  flow.  Figure  1.16  gives  the  pressure  variations 

along  the  liquid- vapor  interface  for  three  different  vapor  flow  models: 

the  present  compressible  and  incompressible  results,  as  well  as  the 

incompressible  similarity  solution  by  Faghri  (1988).  Ve  found  that  the 

results  for  the  incompressible  similarity  solution  and  the  present 

incompressible  numerical  analysis  are  very  close.  The  maximum  deviations 

are  less  than  67,  for  all  three  cases  presented.  The  comparison  between  the 

compressible  and  incompressible  models  shows  a  large  deviation  for  the  case 

of  Re  =  -6.0  (Case  C),  which  is  about  247,  at  the  inlet  of  the  condenser. 

The  maximum  deviations  for  the  case  of  Re  =  -4.0  (Case  B)  and  Re  = 

r  5  c  r  j  c 

-2.0  (Case  A)  are  about  11%  and  37,,  respectively,  with  very  small 
deviations  at  the  end  of  the  condenser.  Therefore,  the  compressibility  of 
the  vapor  does  not  have  a  significant  effect  on  the  total  vapor  pressure 
drop. 

The  outer  wall  temperature  uniformity  is  the  major  concern  for  heat 
pipe  users  and  is  easy  to  measure  experimentally.  The  outer  wall 


67 


sodium  heat  pipe  vapor- liquid 


ect  of  compressibility  on  the  outer  wall  temperature 
variation  along  the  sodium  heat  pipe 


incompressible,  Rer  o=-2.0  compressible.  Re 


variation  along  the  centerline  of  the  sodium  heat  pipe 


temperature  distribution  is  thus  presented  here  in  Fig.  1.17  to  see  the 

effect  of  compressibility.  As  the  results  show,  a  significant  deviation 

exists  between  the  compressible  and  incompressible  models  for  the  case  of 

Re  =  -6.0  (Case  C).  The  incompressible  model  gives  a  more  uniform 

temperature  profile  than  that  of  the  compressible  model.  For  the  case  of 

TQ  -  535°C,  the  largest  deviation  of  these  two  models  is  about  1 0° C .  The 

deviations  are  less  significant  for  the  results  of  Re  =  -4.0  (Case  B) 

and  Re  =  -2.0  (Case  A), 
r ,  e 

Since  the  Mach  number  is  a  measure  of  the  effect  of  the  vapor 
compressibility,  the  variation  of  the  Mach  number,  M,  along  the  centerline 
of  the  pipe  is  shown  in  Fig.  1.18.  The  compressibility  effect  is  important 
when  M  >  0.3  and  is  verified  here  from  the  present  numerical  results  of  the 
heat  pipe  analysis.  Ve  know  from  the  case  of  Re  =  6.0  (Case  C)  that 
when  the  Mach  number  is  less  than  0.2,  very  small  deviations  are  noticed, 
but  when  M  >  0.3,  the  deviation  becomes  significant.  The  maximum  deviation 
of  about  247,  occurs  at  the  exit  of  the  adiabatic  section  with  M=0.66. 
Comparing  Fig.  1.18  with  Figs.  1.16  and  1.17,  we  notice  that  the  deviation 
of  the  models  for  the  pressure  and  temperature  variations  correspond  to  the 
Mach  number  variation.  As  the  Mach  number  increases,  the  deviation  between 
the  compressible  and  incompressible  results  increases. 

Figures  1.19  and  1.20  show  the  effect  of  compressibility  on  the  water 
heat  pipe  vapor  flow  (Case  D) .  A  comparatively  large  radial  Reynolds 
number,  Re  =  -50,  is  chosen  to  see  the  compressibility  effect.  Fig. 
1.19  presents  the  pressure  variation  along  the  vapor- liquid  interface  for 
the  compressible  and  incompressible  models.  The  profiles  almost  overlap 
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compressible  incompressible 


profile  along  the  vapor— liquid  interface  of  the  water  heat  pipe 
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along  the  evaporator  and  condenser  sections  with  only  a  very  small 
deviation  of  less  than  1.57,  in  the  adiabatic  section.  Fig.  1.20  shows  the 
Mach  number  variation  along  the  centerline  of  the  vapor  space.  It  shows 
that  the  maximum  Mach  number  is  only  about  0.018  and  the  results  for  both 
cases  are  very  close  to  each  other.  This  is  because  for  water  as  the 
working  fluid,  the  vapor  density  is  much  greater  than  that  of  sodium  vapor, 
so  the  vapor  velocity  is  much  smaller  for  the  water  heat  pipe  for  a  given 
heat  input.  Actually,  Case  D  has  a  fairly  high  heat  transfer  rate,  but  the 
Mach  number  is  still  very  small.  Thus  we  concluded  that  the  effect  of 
vapor  compressibility  is  not  important  for  heat  pipes  using  water  as  a 
working  fluid. 

1.6.3.  Pressure  Recovery  and  Flow  Reversal 

The  pressure  of  the  vapor  in  the  evaporator  decreases  along  the  path  of 

the  vapor  stream  because  of  friction  and  acceleration  of  the  flow  as  a 

result  of  the  injection  of  vapor.  In  the  condenser  section,  the  extraction 

of  mass  leads  to  a  deceleration  of  the  flow,  i.e.,  to  an  increase  in 

pressure,  while  friction  lowers  the  pressure  of  the  vapor.  The  variation 

of  pressure  in  the  condenser  section  can  be  different,  depending  on  the 

ratio  of  effects  of  friction  and  inertia.  The  results  in  Fig.  1.16  show 

that  the  extent  of  the  pressure  recovery  increases  as  the  condenser  radial 

Reynolds  number  increases.  For  the  case  with  Re  =  1.23  (Case  A),  almost 

no  pressure  recovery  can  be  noticed.  This  is  a  result  of  the  viscous 

effect  dominating  with  the  small  radial  Reynolds  number.  For  Re  =  2.67 

r ,  c 

(Case  B)  and  Re  =4  (Case  C),  the  pressure  recovery  is  approximately  357, 
r ,  c 

and  557,,  respectively.  Fig.  1.19  shows  a  large  pressure  recovery  in  the 
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condenser  section  for  the  heat  pipe  using  water  as  a  working  fluid  (Case 
D).  Up  to  907.  of  the  pressure  drop  in  the  evaporator  and  adiabatic 
sections  were  recovered  at  the  end  of  the  condenser. 

Ve  also  found  that  flow  reversal  takes  place  in  the  condenser  region 

when  the  pressure  recovery  is  significant.  Flow  reversal  was  not  found  for 

the  case  of  Re  =  1.33  (Case  A),  and  a  very  slight  amount  of  flow 

reversal  at  the  end  of  condenser  was  noticed  for  Re  =  2.67  (Case  B), 

which  is  in  agreement  with  what  Bankston  and  Smith  (1973)  observed.  The 

dimensionless  axial  velocity  profile  for  Re  =  4.0  at  three  different 

r>c 

axial  locations  is  plotted  in  Fig.  1.21.  It  shows  clearly  that  the 

velocity  in  the  condenser  section  is  more  extruded  and  flow  reversal  occurs 

near  the  wall  region.  Also,  we  found  for  water  as  the  working  fluid,  a 

larger  Re  compared  with  sodium  as  the  working  fluid  has  to  be  specified 

to  get  flow  reversal.  So  for  different  working  fluids,  the  starting  point 

for  flow  reversal  will  depend  on  different  Re  values. 

r,c 

1.6.4.  Effect  of  Viscous  Dissipation 

The  viscous  dissipation  term  in  the  energy  equation  is  included  in  the 
present  analysis.  Ve  intend  here  to  see  the  effect  of  this  term  on  the 
sodium  heat  pipe  vapor  flow.  For  metallic  working  fluids,  the  density  is 
extremely  small  at  low  vapor  pressures  and  even  for  a  relatively  small  heat 
transfer  rate,  the  vapor  velocity  in  the  axial  direction  can  be  very  large. 
The  results  presented  in  this  section  correspond  to  Case  B  in  Table  1.6  for 
the  compressible  and  incompressible  models.  Fig.  1.22  shows  the  axial 
dimensionless  velocity  profile  in  the  radial  direction  at  three  different 
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compressible  sodium 


The  dimensionless  axial  velocity  profile  at  different  axial  locations 
ncompressible  sodium  heat  pipe  vapor  flow  with  Rere=-4.0  (Case  B) 


locations  for  Case  B  of  the  incompressible  flow.  Only  the  case  with 

viscous  dissipation  is  presented  in  this  figure  since  the  cases  with  and 

without  viscous  dissipation  are  almost  identical.  Since  the  dominant  term 

o 

in  the  viscous  dissipation  ^  of  Eq.  (1.2.4)  is  (dv^/dr)  ,  it  is  helpful  to 
have  a  velocity  profile  for  discussion.  From  Fig.  1.22  the  axial  velocity 
profile  in  the  evaporator  is  very  close  to  that  in  the  adiabatic  section; 
this  means  that  the  velocity  is  nearly  fully  developed  in  the  evaporator 
and  adiabatic  regions.  In  the  condenser  region,  however,  the  mass 
extraction  makes  the  profile  extruded,  which  means  that  it  has  a  larger 
velocity  variation  in  the  radial  direction  which  had  also  been  observed  by 
Busse  (1987a). 

Figure  1.23  shows  the  effect  of  viscous  dissipation  on  the  radial 
temperature  variation  at  three  different  axial  locations  for  incompressible 
vapor  flow  in  the  sodium  heat  pipe  for  Case  B.  The  temperature  variations 
in  the  liquid- wick  and  wall  regions  are  also  presented  in  this  figure. 
Fig.  1.23  also  shows  that  the  viscous  dissipation  can  change  the  radial 
temperature  distribution  significantly.  In  the  beginning  of  the 

evaporator,  the  two  temperature  curves  are  quite  close,  but  in  the 
adiabatic  and  condenser  regions,  the  temperature  difference  .  rt/. een  the 
curves  which  include  and  do  not  include  viscous  dissipatio*  are  much 
larger.  For  this  particular  case  with  Tq  =  535°C,  it  can  be  as  much  as  8°C 
difference.  Also,  the  viscous  dissipation  term  increases  the  temperature 
more  near  the  wall  because  of  the  larger  velocity  gradient,  as  shown  in 
Fig.  1.22. 
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Figure  1.24  presents  the  effect  of  viscous  dissipation  on  the  axial 
temperature  variation.  In  PHOENICS,  only  the  viscous  dissipation  term  for 
incompressible  flow  can  be  controlled  using  one  command.  This  is  because 
another  term  DP/Dt  is  also  involved  in  the  same  command  for  t.he 
compressible  flow  model.  Ve  cannot  deactivate  the  viscous  term  alone. 
Therefore,  the  results  are  discussed  based  on  the  incompressible  flow 
model,  but  the  results  with  viscous  dissipation  for  the  compressible  model 
are  also  presented  for  the  axial  temperature  variations.  The  figure  shows 
that  for  incompressible  flow,  the  interface  temperature  almost  remains  the 
same  both  with  and  without  the  viscous  dissipation  terms.  However,  a  large 
difference  exists  in  the  vapor  temperature  along  the  centerline  of  the 
pipe.  In  the  evaporator,  the  two  cases  are  close  to  each  other  but  then 
start  to  deviate  significantly  from  each  other  in  the  adiabatic  section. 
This  deviation  continues  through  the  heat  pipe  to  the  end  of  the  condenser 
section.  The  largest  deviation  occurs  near  the  end  of  the  condenser.  The 
centerline  temperature  for  compressible  flow  is  higher  than  that  of  the 
incompressible  model,  which  means  that  the  viscous  term  has  a  more 
significant  effect  on  the  compressible  flow  model.  This  is  probably 
because  the  momentum  equation  is  coupled  with  the  energy  equation  by  the 
density  variation,  so  the  viscous  term  will  affect  the  pressure  and 
velocity  profiles  and  thus  the  overall  temperature  distribution.  For  the 
incompressible  model,  the  only  link  between  the  momentum  and  energy 
equations  is  the  vapor  heat  conduction  term  at  the  vapor- liquid  interface 
in  Eq.  (1.15),  but  the  effect  is  very  small.  Therefore,  the  velocity  fields 
with  and  without  the  viscous  terms  are  almost  identical  to  each  other. 
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viscous  dissipation  on  the  axial  temperature 
sodium  heat  pipe  with  Rer  ~=-4.0  (Case  B) 


Since  the  temperature  boundary  condition  at  the  interface  is  controlled  by 
the  Clapeyron  correlation,  the  interface  temperature  variations  are 
slightly  different  for  both  cases. 

1.6.5.  Elliptic  and  Parabolic  Comparison 

Figure  1.25  presents  a  comparison  between  the  normalized  pressure  drop 

profiles  of  the  elliptic  and  partially  parabolic  solutions  for  the 

compressible  heat  pipe  vapor  flow.  For  the  partially  parabolic  solution, 

the  axial  diffusion  terms  are  neglected  in  the  governing  equations.  The 

results  show  that  the  partially  parabolic  solutions  have  larger  pressure 

drops  and  the  deviation  ranges  from  1-117.  for  the  cases  presented.  For  the 

case  using  water  as  the  working  fluid  with  a  higher  radial  Reynolds  number 

of  Re  =  -50  (Case  D),  a  larger  deviation  of  117.  is  noticed  at  the  end  of 

the  condenser.  For  the  case  of  a  sodium  heat  pipe  with  a  radial  Reynolds 

number  of  Re  -  -6.0  (Case  C),  the  deviation  is  about  77,.  Vhen  the 

r,e 

radial  Reynolds  number  is  decreased  to  Re  =  -2.0,  the  two  curves  almost 

r 

overlap.  This  trend  is  similar  to  what  Tien  and  Rohani  (1974)  observed, 
but  the  differences  found  here  between  these  two  models  are  not  as  large  as 
what  they  observed.  Faghri  and  Parvani  (1988)  did  a  comparison  between  the 
elliptic  and  partially  parabolic  solutions  for  the  incompressible  vapor 
flow  in  a  concentric  annular  heat  pipe.  In  that  study,  very  small 
deviations  were  observed  due  to  the  incompressible  flow  analysis.  However, 
the  present  results  show  that  it  is  better  to  use  the  elliptic  approach 
when  the  radial  Reynolds  number  is  large  for  the  compressible  model. 


82 


-50 


iparison  between  the  elliptic  and  parabolic  solutions  for  the  axial 
sure  profile  along  the  heat  pipe  vapor- liquid  interface 


Since  no  axial  heat  conduction  through  the  wall  is  considered  for  the 
present  partially  parabolic  solution,  the  total  heat  input  at  the 
liquid- vapor  interface  for  the  parabolic  solution  is  greater  than  that  of 
the  elliptic  solution.  As  a  result,  more  liquid  is  evaporated  and  a  larger 
pressure  drop  is  expected  for  the  partially  parabolic  solution.  As 
expected,  for  the  case  with  a  high  discussed  in  the  section  concerning 
the  conjugate  effect,  the  elliptic  and  partially  parabolic  solutions  have 
larger  deviations  because  of  the  significant  axial  heat  conduction  through 
the  pipe  wall.  The  maximum  deviation  of  the  pressure  variation  between  the 
two  solutions  for  =  1000  (Case  B3)  is  about  127.. 
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1.7  Conclusions 


The  high  and  low  temperature  heat  pipes  have  been  studied  numerically  by 
using  the  generalized  computer  code  PHGENICS  for  steady- state  operation. 
The  objective  of  this  study  is  to  find  the  effects  of  conjugate  heat 
transfer  and  vapor  compressibility.  The  following  conclusions  have  been 
made: 

1.  The  numerical  results  show  a  fairly  good  agreement  compared  with  the 

experimental  data  at  both  low  and  high  operating  temperatures.  It  is 
also  believed  that  the  present  model  can  predict  the  general 

performance  of  heat  pipes  with  single  or  multiple  heat  sources. 

2.  The  axial  conduction  through  the  wall  and  liquid- wick  with  large  values 

of  has  a  positive  effect  on  the  heat  pipe  performance.  The 

conjugate  axial  conduction  tends  to  make  the  outer  wall  and  interface 
temperature  more  uniform.  For  the  case  with  relatively  small  values  of 
Kwj  and  A,  the  results  approach  the  solution  of  the  case  when  there  is 
no  axial  conduction  considered. 

3.  (a)  The  vapor  compressibility  can  be  important  for  the  prediction  of 
the  sodium  heat  pipe  temperature  profile  when  the  Mach  number  is 
greater  than  0.3. 

(b)  The  compressible  and  incompressible  models  predict  almost  the  same 
value  for  the  total  pressure  drops. 
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(c)  For  the  incompressible  model,  the  similarity  solution  by  Faghri 
(1988)  shows  a  very  good  agreement  (less  than  1.57.  deviation)  for  the 
pressure  profile  prediction  compared  with  the  present  numerical 
results. 

4.  Up  to  907.  of  the  pressure  recovery  in  the  condenser  section  has  been 

found  for  the  heat  pipe  using  water  as  a  working  fluid.  The  flow 

reversal  for  a  sodium  heat  pipe  vapor  flow  begins  at  Re  =  2.67  (Case 

r ,  c 

B)  for  the  cases  studied. 

5.  The  viscous  dissipation  term  in  the  energy  equation  does  affect  the 
temperature  profile  inside  the  vapor  space,  but  has  a  very  small  effect 
on  the  interface  temperature  profile. 

6.  In  general,  the  partially  parabolic  solutions  are  very  close  to  the 
elliptic  solutions  except  for  the  cases  with  a  large  radial  Reynolds 
number  at  the  interface  or  large  values  of  K^. 

7.  For  the  heat  pipe  using  water  as  the  working  fluid,  the  vapor 
temperature  along  the  heat  pipe  is  almost  uniform.  This  is  due  to  a 
very  small  pressure  drop  along  the  pipe  compared  with  the  static  vapor 
pressure  and  also  the  thermodynamic  equilibrium  between  the  pressure 
and  temperature  at  the  liquid- vapor  interface. 
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II.  SIMULTANEOUS  AXIAL  CONDUCTION  IN  THE  FLUID  AND  THE  PIPE 
FOR  FORCED  CONYECTIVE  LAMINAR  FLOV  VITH  BLOVING 
AND  SUCTION  AT  THE  WALL 

2 . 1  Summary 

Numerical  solutions  are  reported  for  conjugate  heat  transfer  in  a 
porous  pipe  having  an  internal  laminar  flow  with  blowing  or  suction  at 
the  inner  surface  of  the  pipe  and  constant  heat  flux  at  the  outer  surface 
of  the  pipe.  The  effect  of  the  simultaneous  axial  conduction  through 
the  wall  and  the  fluid  has  been  studied  for  the  combined  hydrodynamic  and 
thermal  entry  lengths.  The  results  show  that  the  ratio  of  the  thermal 
conductivities  of  the  pipe  wall  to  the  fluid  and  the  thickness  of  the 
pipe  wall  may  become  significant  factors  on  the  heat  transfer  when  the 
Peclet  number  is  small,  especially  for  the  case  when  fluid  is  injected 
into  the  pipe.  Also,  the  effect  of  axial  wall  conduction  for  the  case  of 
constant  heat  flux  at  the  outer  wall  surface  can  be  neglected  when  the 
wall  thickness  is  small  and  the  ratio  of  the  conductivities  of  the  wall 
to  the  fluid  approaches  unity. 
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2.2  Introduction 


Conjugate  heat  transfer  in  a  pipe  with  internal  laminar  fluid  flow 
is  a  problem  of  considering  the  simultaneous  heat  transfer  inside  the 
fluid  and  the  solid  wall.  In  conventional  convective  heat  transfer 
problems,  the  thermal  boundary  condition  at  the  solid- fluid  interface  is 
assumed  to  be  known,  either  in  terms  of  the  heat  flux  or  the  temperature. 
Some  considerations  have  been  given  in  the  past  concerning  the  errors 
that  may  be  introduced  by  this  assumption. 

The  primary  research  on  the  conjugate  heat  transfer  problem  for  the 
circular  tube  was  carried  out  by  Luikov  et  al.  (1971),  but  their 
closed- form  solution  involved  highly  complicated  functions.  Because  of 
this  fact,  no  numerical  results  were  reported.  A  detailed  solution  was 
given  by  Mori  et  al.  (1974,  1976)  for  constant  heat  flux  and  constant 
temperature  at  the  outer  surface  of  the  pipe.  The  solution  was  obtained 
by  assuming  that  the  velocity  profile  was  the  fully  developed  parabolic 
profile  and  by  neglecting  the  axial  conduction  of  the  fluid.  Mori  et  al. 
(1974)  also  assumed  the  wall- fluid  interfacial  temperature  distribution 
in  the  axial  direction  was  in  the  form  of  a  power  series  with  unknown 
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coefficients.  With  this  temperature  distribution  as  the  interfacial 
boundary  condition,  the  analytical  solution  to  the  energy  equation  for 
the  fluid  was  obtained  using  the  Graetz  solutions.  The  solution  to  the 
energy  equation  for  the  wall  was  then  derived  readily  since  all  the 
boundary  conditions  were  now  known.  Barozzi  and  Pagliarini  (1985) 
employed  a  finite  element  method  in  conjuction  with  Duhamel’s  theorem  to 
extend  the  results  given  by  Mori  et  al.  (1974)  for  different  pipe 
lengths.  Results  by  Mori  (1974),  and  Barozzi  and  Pagliarini  (1985)  are 
applicable  to  cases  with  short  heating  sections  due  to  the  insulated 
boundary  conditions  imposed  on  the  wall  upstream  and  downstream  of  the 
heated  section. 

The  coupled  effect  of  axial  conduction  in  the  wall  and  the  fluid  in 
laminar  pipe  flow  with  an  applied  wall  heat  flux  in  the  downstream  region 
was  studied  for  very  long  ducts  numerically  by  Faghri  and  Sparrow  (1980) 
and  Zariffeh  et  al.  (1982)  and  analytically  by  Campo  and  Rangel  (1983) 
and  Soliman  (1984).  In  these  analyses  the  pipe  is  extended  indefinitely 
upstream  and  downstream  of  the  heating  section.  Faghri  and  Sparrow 
(1980)  investigated  simultaneous  wall  and  fluid  axial  conduction  in 
laminar  pipe  flow  with  the  assumption  that  the  pipe  wall  was  thin  and  the 
velocity  profile  was  the  fully  developed  parabolic  profile.  They 
observed  that  axial  conduction  inside  the  tube  wall  can  cause  a 
substantial  preheating  of  both  the  wall  and  fluid  in  the  upstream  region 
where  there  is  no  external  heat  input.  Conjugate  heat  transfer  when  the 
fluid  is  turbulent  was  investigated  by  Kuznetzov  and  Belousov  (1974), 
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Sakakibara  and  Endoh  (1977)  and  Lin  and  Chow  (1984). 

Recently,  more  attention  has  been  given  to  the  fluid  flow  and  heat 
transfer  in  porous  pipes  with  applications  such  as  transpiration  cooling 
and  heat  pipes.  Usually,  the  porous  pipe  wall  is  of  a  significant 
thickness,  so  that  axial  wall  conduction  may  have  an  influence  on  the 
heat  transfer.  In  a  heat  pipe,  the  thin  wick  is  attached  to  the  inner 
wall  of  a  thick  tube  to  achieve  the  capillary  force  for  the  liquid 
return.  Blowing  and  suction  occurs  because  of  the  evaporation  and 
condensation  of  the  working  fluid  in  the  heating  and  cooling  segments  of 
the  heat  pipe.  This  application  was  the  motivation  for  the  present  work. 
No  literature  has  been  found  concerning  the  problem  of  conjugate  heat 
transfer  with  blowing  and  suction  at  the  wall.  The  objective  of  the 
present  analysis  is  to  investigate  the  effect  of  axial  conduction  in 
porous  pipes  in  the  combined  hydrodynamic  and  thermal  entry  lengths  with 
constant  heat  flux  at  the  outer  surface  of  the  pipe  wall.  In  addition, 
the  influence  of  the  Peclet  number  is  stressed  due  to  the  inclusion  of 
the  axial  conduction  term  in  the  fluid  so  as  to  complete  the  analysis 
given  by  Mori  et  al.  (1974).  The  numerical  solution  relies  on  the 
control  volume  formulation  by  Patankar(1980)  and  an  iterative  scheme 
which  dealt  simultaneously  with  the  fluid  and  the  wall.  The  numerical 
results  for  the  Nusselt  number  and  temperature  at  the  interface  are 
presented  as  in  by  Mori  (1974,  1976).  In  addition,  the  interfacial  heat 
flux  for  both  porous  and  impermeable  wall  cases  are  also  given  to  make 
the  results  useful  in  engineering  applications. 
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2.3  Mathematical  Formula! 4 "n 

The  model  used  for  analysis  is  shown  in  Figure  2.1.  A  constant  heat 
flux  qQ  is  applied  uniformly  along  the  outer  surface  of  the  wall  over  a 
finite  length  L.  The  fluid  enters  the  inlet  of  the  pipe  with  a  uniform 
velocity  w^n  and  uniform  temperature  T^.  The  injected  and  extracted 
fluid  is  the  same  as  that  of  the  main  pipe  flow  and  is  at  the  interface 
temperature  when  it  enters  or  leaves  the  walls.  The  equations  governing 
the  present  problem  are  the  conservation  of  mass,  momentum  and  energy 
equations  whereby  assuming  laminar,  steady  and  incompressible  flow  with 
constant  fluid  properties  are  reduced  to  the  following  forms  in 
cylindrical  coordinates. 


(2.1) 


(2.2) 


(2.3) 


(2.4) 
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■transfer  section  of  a  porous  pipe 


Ve  should  note  that  in  equations  (2.2- 2.4),  the  terms  in  braces  {} 
are  associated  with  axial  diffusion  terms.  These  terms  are  neglected 
when  the  partially  parabolic  version  is  considered  but  are  accounted  for 
in  the  elliptic  case.  The  thermal  conductivity  k  of  the  solid  is  in 
general  different  from  that  of  the  fluid.  At  the  solid- fluid  interface, 
the  harmonic  mean  is  used  to  determine  the  thermal  conductivity. 

The  boundary  conditions  for  the  case  of  uniform  inlet  conditions  and 
constant  heat  flux  at  the  outer  wall  are  defined  as  follows: 


(2.5) 

(2.6) 

(2.7) 

(2.8) 
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v^  q  suction 

v  =  'i 

v^  <  0  blowing 

(2.9) 

T,  =  T 
f  w 

v^  =  0  impermeable  wall 

(2.10) 

0Tf 

*rw 

kf  W  = 

kw  w 

(2.11) 

Outlet  of  the  tube  (elliptic  solution):  z  =  L 

91 

ri  <  r  <  V  W~  =  0  (2-12) 

6Tr 

0  <  t  <  t i?  =  0,  p  =  0  (2.13) 

Ve  should  emphasize  that  the  physical  situation  corresponding  to  this 
problem  is  the  one  in  which  conduction  is  confined  to  a  finite  length  of 
heated  wall,  with  an  adiabatic  boundary  condition  imposed  at  the  upstream 
and  downstream  ends  of  the  heated  segment.  This  constraint  restricts  the 
results  to  tubes  with  short  heating  sections. 

The  governing  equations  (2. 1-2.4)  with  the  appropriate  boundary 
conditions  (2.5-2.13)  contain  five  independent  parameters: 

*  Peclet  number  of  the  fluid, 

Pe  =  Re  *  Pr  =  2  r^  win  pf  Cp/kf 
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*  Radial  Reynolds  number  at  the  solid- fluid  interface, 


*  Vail- to- fluid  thermal  conductivity  ratio, 

K  -  k„/kf 

*  Thickness  of  the  wall, 

A  =  (Vri)/2ri 

*  Heated  length, 

L*  =  L/2riPe 


The  computations  are  carried  out  for  the  Peclet  number  of  the  fluid 

of  100  and  1000  with  A  of  0.01  and  0.1  and  K  =  1,  500  and  5000  for 

different  suction  and  blowing  radial  Reynolds  numbers  at  the  wall.  This 

choice  of  parameters  covers  a  wide  range  of  possible  combinations  of 

fluid  and  pipe  vail  properties,  flow  rates  and  boundary  condition 

specifications.  The  negative  and  positive  radial  Reynolds  numbers  at  the 

wall  indicate  blowing  and  suction,  respectively. 

The  numerical  results  can  be  presented  for  the  following  four 

* 

dimensionless  groups  in  terms  of  the  dimensionless  axial  distance  z  = 
z 

"SrTPe  : 


Interface  temperature, 


T,  -  T. 


9 i ■ 


in 

~E7rf 
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Bulk  Temperature, 
T.  -  T- 

a  _  b  in 


%  vkf 


*  Nusselt  number, 

Nuz =  T-- t: 

i  b 


where 


(M)  -  ** 

Wr=r.“  "ij 

*  Heat  flux  across  the  interface, 
^i/^mean 


Although  the  Nusselt  number  is  traditionally  the  main  dimensionless 
parameter  used  in  presenting  results  for  conventional  convection 
problems,  there  is  a  good  justification  not  to  do  so  for  conjugate  heat 
transfer  problems  because  q^,  T^,  are  all  unknowns  in  the  definition 
of  the  Nusselt  number  as  given  above.  Mori  et  al.  (1974)  presented 
results  only  for  Nu  and  9 ■ ,  and  therefore  their  results  are  of  no  use 
for  direct  engineering  applications.  From  the  engineering  viewpoint,  the 
most  important  information  is  the  rate  of  heat  transfer  across  the 
interface.  So  the  results  are  presented  in  terms  of  the  dimensionless 
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group  q./q  as  well  as  the  Nusselt  number.  The  local  heat  flux  across 
o  r  'iniean 

kh  (T.  -  T.  f) 

the  interface  q^  is  approximated  by  'r',in  j\\ — y  with  kh  being  the 

harmonic  mean  of  the  thermal  conductivities  of  the  solid  wall  and  the 
fluid  which  is  defined  as  2k£-kw/(k£+kw)  by  Patankar  (1978).  The  radial 
distance  from  the  pipe  axis  to  the  interface  is  r-  and  T.  and  T.  ,  are 

1  1  j  n  i.  j  1 

the  temperatures  at  the  grid  nodes  adjacent  to  the  interface  on  the  wall 

and  the  fluid  sides,  respectively.  The  constant  mean  heat  flux  at  the 

interface  without  considering  axial  heat  conduction  in  the  tube  wall 

q  is  equal  to: 
unean  n 


qmean  ~  %  r^ 
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2.4  Solution  Procedure 


The  problem  is  solved  as  a  convection- conduction  problem  throughout 
the  entire  calculation  domain,  but  since  the  velocities  in  the  solid  are 
zero,  the  analysis  in  the  solid  region  was  for  a  pure  conduction  problem. 
To  account  for  the  discontinuity  in  the  value  of  the  thermal  conductivity 
at  the  wall- fluid  interface,  the  transport  coefficient  in  the  energy 
equation  is  evaluated  as  the  harmonic  mean  of  the  values  on  each  side  of 
the  interface  developed  by  Patankar  (1980). 

The  finite  difference  iteration  method  of  solution  developed  by 
Spalding  et  al.  (1980)  which  is  employed  in  the  generalized  PHOENICS 
Computational  Code  by  Spalding  and  Rosten  (1985)  is  used  for  the  elliptic 
and  partially  parabolic  solutions  of  equations  (2. 1-2.4)  with  the 
boundary  specifications  given  in  equations  (2.5-2.13).  By  this  method, 
the  above  equations  are  solved  over  a  square  mesh  by  the  finite  control 
volume  method  outlined  by  Patankar  (1980).  Finite  control  volume 
equations  are  derived  by  the  integration  of  the  differential  equations 
over  an  elementary  control  volume  or  a  cell  surrounding  a  grid  node. 
Upwind  differencing  is  used  for  the  convective  terms. 

Pressures  are  solved  from  a  pressure  correction  equation  which 
yields  the  pressure  change  needed  to  procure  velocity  changes  to  satisfy 
mass  continuity.  The  ’SIMPLEST’  practice  by  Spalding  (1980)  is  employed 
for  the  momentum  equations.  The  most  significant  difference  between 
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’SIMPLEST’  and  the  ’SIMPLE’  algorithm  by  Patankar  (1980)  is  that  in  the 
former  the  control  volume  formulation  coefficients  for  momentum  contain 
only  diffusion  contributions,  the  convection  terms  being  added  to  the 
linearized  source  term  of  the  equations. 

The  equations  are  solved  by  a  line-by-line  procedure  which  is 
similar  to  Stone’s  Strongly  Implicit  Method  but  is  free  from  parameters 
requiring  case- to- case  adjustment  and  is  therefore  less  complex.  The 
temperature  is  solved  in  a  whole- field  manner,  and  the  pressure 
correction  equation  is  solved  in  a  slab- by- slab  manner. 

The  accuracy  of  the  numerical  solution  was  checked  in  two  ways:  the 
grid  spacing  was  systematically  varied  and  the  results  for  different  grid 
sizes  were  compared  to  the  extrapolated  results  of  the  infinitesimal  grid 
spacing,  and  the  numerical  results  were  compared  with  a  known  solution 
for  uniform  heat  flux  and  temperature  at  the  interface  without  axial 
conduction  in  the  wall  or  the  fluid  and  with  no  blowing  or  suction  by 
Hornbeck  (1966).  The  test  was  passed  at  a  satisfactory  level  of 
agreement,  within  17,  for  most  of  the  results. 

A  number  of  different  uniform  grid  sizes  were  chosen  to  test  the 
accuracy  of  the  solution.  After  considering  the  accuracy,  computer  time 
and  the  truncation  error,  the  final  grid  sizes  that  were  chosen  for  the 
presentation  of  results  are  as  follows: 

20x4x40  (radial  fluid  x  radial  solid  x  axial)  for  A  =  0.01 

20x10x20  (radial  fluid  x  radial  solid  x  axial)  for  A  =  0.1 


99 


A  converged  solution  was  obtained  by  checking  the  variation  of  the 
residuals  in  terms  of  the  sweep  number.  When  increasing  the  sweep  number 
to  the  extent  that  no  further  decrease  of  the  residuals  is  observed,  this 
means  the  results  are  converging.  For  this  conjugate  problem,  150  sweeps 
were  adequate  for  convergence  of  the  elliptic  solution. 
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2.5  Results  and  Discussion 

The  numerical  results  were  obtained  for  at  least  two  different 

values  of  each  of  the  dimensionless  parameters  in  order  to  determine  the 

importance  of  axial  conduction.  These  parameters  are  the  Peclet  number 

of  the  fluid,  the  radial  Reynolds  number  at  the  solid- fluid  interface, 

the  wall- to- fluid  thermal  conductivity  ratio,  and  the  dimensionless 

thickness  of  the  pipe  wall.  In  the  present  work,  only  one  dimensionless 
* 

pipe  length  L  =  0.05  is  chosen  to  avoid  the  excessive  presentation  of 

the  results.  Axial  heat  conduction  along  the  pipe  wall  is  usually 

neglected  by  practicing  engineers  when  designing  heat  exchangers  or  other 

kinds  of  heat  transfer  devices  because  it  is  very  difficult  to  deal  with 

conjugate  heat  transfer.  Therefore,  it  is  of  practical  importance  to 

know  the  errors  that  may  be  introduced  by  this  assumption.  The  numerical 

results  for  the  three  dimensionless  variables:  the  Nusselt  number  Nu  , 

z 

the  interface  temperature  9^,  and  the  heat  flux  across  the  interface 

* 

q./q  are  given  in  terms  of  z  for  different  values  of  K,  A,  Pe  and 
Hr  unean  b  ’ 

Re^.  The  bulk  temperature  0^  can  be  calculated  from  the  information 
presented. 

Figures  2.2  through  2.4  show  the  effect  of  the  thermal  conductivity 
ratio  for  A  =  0.1,  Pe  =  100  and  Re^  =  1.0,  0.0,  -3.0  on  the  axial 
variation  of  the  Nusselt  number,  the  interface  temperature  and  the 
interfacial  heat  flux,  respectively.  The  effect  of  axial  conduction 
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Fig.  2.2(b)  The  effect  of  the  thermal  conductivity  ratio  K  on  the 
Nusselt  number  at  the  interface  for  Pe=100  and  A=0.1 
with  suction  and  blowing 


u 


^rab//!b 
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Fig.  2.4  The  effect  of  the  thermal 
the  interfacial  heat  flux  for 


becomes  more  important  as  the  conductivity  ratio,  K,  increases.  As  K 

increases,  the  results  show  that  the  interface  temperature  distribution 

* 

changes  gradually  to  a  constant  temperature  along  z  ,  and  the  Nusselt 
number  approaches  that  of  a  constant  wall  temperature.  Physically,  when 
K  is  very  large  the  axial  conduction  in  the  wall  has  a  tendency  to  make 
the  wall  temperature  uniform.  It  is  similar  to  the  situation  when 
condensation  or  evaporation  at  a  wall  occurs,  which  leads  to  a  uniform 
wall  temperature  because  of  the  high  heat  transfer  coefficient.  The 
blowing  case  has  a  higher  interface  temperature  and  the  suction  case  has 
a  lower  interfacial  temperature  in  comparison  with  the  case  in  which  the 
wall  is  impermeable.  At  K  =  5000,  the  interface  temperature  distribution 
becomes  constant  except  for  the  short  entrance  region  and  the  Nusselt 
number  also  approaches  the  constant  temperature  results  by  Hornbeck 

(1966).  It  is  noticed  from  the  present  results  and  those  presented  by 
Raithby  (1971)  for  the  case  of  no  axial  conduction  that  blowing  at  the 
wall  caused  a  bigger  difference  between  the  Nusselt  numbers  for  the  cases 
of  constant  heat  flux  and  constant  wall  temperature.  Sr>  the  effect  of 
axial  conduction  is  more  pronounced  in  this  case  than  in  other  cases.  In 
general,  as  the  thermal  conductivity  decreases,  the  Nusselt  number 

increases  and  the  interfacial  heat  flux  tends  to  become  more  uniform. 

For  K  =  1,  the  results  show  that  the  interfacial  heat  flux  almost  becomes 
uniform  for  all  three  cases.  Also,  the  Nusselt  number  approaches  the 

case  of  constant  heat  flux  at  the  interface.  This  is  because  the  axial 
conduction  causes  both  the  interfacial  temperature  and  bulk  temperature 
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to  increase  at  approximately  the  same  rate,  so  that  the  difference  of 
T^-T^  and  ($j)r_r>  almost  remains  unchanged.  The  numerical  results  for 

constant  heat  flux  and  temperature  at  the  interface  for  blowing  and 
suction  are  computed  with  the  consideration  of  axial  conduction  through 
the  fluid.  The  numerical  results  for  the  Nusselt  number  for  K  =  1 
without  blowing  and  suction  compares  favorably  with  the  numerical  results 
for  the  case  of  constant  heat  flux  at  the  interface  as  given  by  Hornbeck 
(1966).  The  Nusselt  number  for  K  =  1  with  suction  is  slightly  smaller 
than  the  results  for  constant  heat  flux  at  the  interface  and  the  results 
for  the  case  of  blowing  is  slightly  greater  than  the  results  for  constant 
heat  flux  at  the  interface.  This  is  probably  because  the  results  for 
constant  heat  flux  at  the  interface  include  the  axial  heat  conduction 
through  the  fluid.  The  results  for  the  Nusselt  number  with  K  >  1  are 
between  two  limiting  cases:  constant  wall  temperature  and  constant  heat 
flux  at  the  interface.  This  behavior  was  first  noticed  by  Mori  et  al. 
(1974)  for  impermeable  walls.  The  present  analysis  shows  also  the  same 
general  trend  for  the  cases  of  blowing  and  suction.  It  appears  that  the 
influence  of  axial  wall  heat  conduction  must  be  accounted  for  if  the 
value  of  K  is  large. 

Figures  2.5  through  2.7  show  the  effect  of  varying  the  tube  wall 
thickness  on  the  Nusselt  number,  the  interface  temperature  and  the 
interfacial  heat  flux,  respectively.  From  these  results,  we  found  that 
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const,  heat  flux  at  interface  (Hornbeck,  1966) 
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Fig.  2.5(a)  The  effect  of  wall  thickness  A  on  the  Nusselt 
number  at  the  interface  for  Pe=100,  K=500  and  RepO.O 


i(b)  The  effect  of  wall  thickness  A  on  the  Nusselt  number 
interface  for  Pe=100  and  K=500  with  suction  and  blowing 


for  small  A  the  effect  of  axial  heat  conduction  through  the  wall  can  be 
neglected  since  the  results  approach  that  of  constant  heat  flux  at  the 
interface.  Faghri  and  Sparrow  (1980)  assumed  a  thin- walled  tube,  so  that 
all  of  the  results  for  fully  developed  flow  approach  the  results  for  a 
constant  heat  flux  at  the  interface.  Ve  found  that  when  the  wall 
thickness  increases,  the  axial  wall  conduction  becomes  more  important, 
and  the  interface  temperature  becomes  more  uniform  because  axial  heat 
conduction  tends  to  level  out  the  temperature  along  the  pipe.  These 
results  also  show  that  the  interfacial  heat  flux  tends  to  be  uniform  when 
the  wall  becomes  thinner.  It  is  interesting  to  notice  that  all  of  the 
results  for  different  values  of  A  are  within  the  two  limiting  cases: 
constant  heat  flux  at  the  interface  and  constant  wall  temperature  at  the 
interface.  Ve  should  mention  that  when  K  1,  the  axial  conduction 
becomes  less  important  so  that  the  influence  of  A  will  also  become  less 
important.  Vhen  K  -*  od  the  influence  of  A  will  become  pronounced.  This 
seems  quite  reasonable  from  the  physical  viewpoint  of  the  problem. 
Apparently  from  these  results,  we  should  pay  more  attention  to  the 
blowing  case  with  large  K.  The  trends  of  the  results  concerning  the 
influence  of  the  wall  thickness  agrees  fairly  well  with  those  by  Mori 
(1974)  for  the  case  of  an  impermeable  wall. 

Figures  2.8  through  2.10  show  the  comparison  between  the  partially 
parabolic  and  the  elliptic  solution  on  the  Nusselt  number  for  different 
values  of  Pe  and  K.  Figure  2.8  shows  the  values  of  the  Nusselt  number 


112 


113 


Fig.  2.8  The  effect  of  axial  heat  conduction  on  the  Nusselt 
number  at  the  interface  for  Pe=100,  K=500  and  A=0.1 
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Fig.  2.10  The  effect  of  axial  heat  conduction  on  the  Nusselt 
number  at  the  interface  for  Pe=1000,  K=500  and  A=0.1 


predicted  with  and  without  wall  and  fluid  axial  conduction  for  Pe  -  100, 
A  =  0.1,  K  =  500  and  Re^  =  1.0,  0.0,  -3.0.  When  axial  conduction  is 
neglected,  the  Nusselt  number  is  overpredicted  along  the  entire  entry 
region.  This  effect  diminishes  as  the  wall- to- fluid  conductivity  ratio 
is  reduced  as  shown  in  Figure  2.9.  When  the  Peclet  number  is  increased 
to  1000,  smaller  differences  are  observed. 
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2.6  Conclusions 


The  effect  of  axial  heat  conduction  have  been  studied  for  cases  with 
blowing  and  suction  and  with  an  impermeable  wall.  Ve  concluded  that 
large  errors  will  arise  in  the  solution  when  the  axial  wall  conduction  is 
neglected  for  blowing  cases  with  large  K.  As  K  increases,  the  results 
approach  the  solution  of  constant  wall  temperature  at  the  interface 
without  axial  wall  conduction.  For  K  =  1,  the  local  Nusselt  number 
approaches  the  solution  of  constant  heat  flux  at  the  interface.  It  is 
confirmed  that  the  effect  of  axial  conduction  in  the  wall  can  be 
neglected  reasonably  when  the  wall  is  very  thin.  But  when  the  wall  is 
very  thick,  the  results  approach  the  solution  of  constant  wall 
temperature  at  the  interface.  Increasing  Pe  by  a  factor  of  10  reduces 
the  effect  of  axial  conduction  much  more  than  a  corresponding  reduction 
of  A  or  K. 
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III.  THE  THERMAL  PERFORMANCE  OF  HEAT  PIPES  VITH  LOCALIZED  HEAT  INPUT 

3.1  Summary 

The  performance  of  heat  pipes  with  localized  heat  input  including 
the  effects  of  axial  and  circumferential  heat  conduction  under  high  and 
low  working  temperatures  was  investigated.  The  numerical  results  show 
that  when  heat  pipes  are  spot  heated,  the  peak  temperature  of  the  wall  is 
greatly  reduced  and  the  surface  can  be  protected  from  being  burned  out  by 
the  high  heat  flux.  The  boiling  limitation  becomes  the  most  important 
limitation  for  this  type  of  heat  pipe.  Numerical  results  for  block 
heating  a  heat  pipe  with  low  working  temperatures  indicate  a  good 
agreement  with  existing  experimental  data.  Also,  most  of  the  input  heat 
passes  through  the  wall  beneath  the  heated  block. 


118 


3.2  Introduction 


Since  the  publication  of  the  first  paper  on  heat  pipes,  various 
kinds  of  heat  pipes  have  been  manufactured,  tested,  and  put  into 
operation.  In  the  meantime,  thousands  of  theoretical  and  experimental 
analyses  dealing  with  the  characteristics  of  heat  pipes  have  been 
published.  Among  the  special  types  of  heat  pipes,  localized  heat  input 
or  spot  heated  heat  pipes  have  not  been  extensively  studied.  This  is 
suprising  since  many  high  performance  heat  pipes  are  subjected  to 
localized  heating  for  a  variety  of  applications.  The  study  of  spot 
heating  heat  pipes  is  important  in  research  areas  such  as  the  cooling  of 
leading  edges  on  hypersonic  aircraft,  the  protection  of  special  surfaces 
from  being  attacked  by  very  high  heat  flux  sources  such  as  a  laser  beam, 
cooling  of  microelectronic  elements,  etc. 

According  to  the  working  conditions  and  the  application,  spot  heated 
heat  pipes  can  be  classified  into  two  major  types:  namely,  heat  pipes 
with  low  or  moderate  working  temperatures  mainly  used  for  the  purpose  of 
energy  conservation  or  electronic  cooling,  and  heat  pipes  with  high 
working  temperatures  used  to  protect  a  surf  xe  from  being  burned  out  by  a 
very  high  heat  flux.  Rosenfeld  [1987]  studied  the  performance  of  line 
heated  heat  pipes  with  low  working  temperatures  analytically  with  a 
one-dimensional  model  (circumferential  direction)  and  numerically  with  a 
two-dimensional  model  (radial  and  circumferential  directions). 

For  the  present  analysis  of  spot  heating  or  block  heating,  we  should 
consider  axial  conduction,  which  has  a  much  more  pronounced  effect  than 
conduction  in  the  radial  direction.  Furthermore,  the  effect  of  radiation 
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is  included  in  the  analysis  of  spot  heated  heat  pipes  with  high  working 
temperatures.  In  addition,  the  operating  temperature  of  the  heat  pipe 
should  be  obtained  by  an  overall  energy  balance  rather  than  an  input 
condition  as  done  by  Rosenfeld  [1987]  for  the  line  heated  heat  pipe. 
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3.3  Analysis 


3.3.1  Spot  heated  heat  pipes  with  high  working  temperatures 

Spot  heated  heat  pipes  should  be  used  to  protect  a  surface  from 
being  burned  out  by  a  very  high  heat  flux  by  using  the  working  fluid 
within  the  porous  media  inside  the  pipe  as  an  evaporator  to  absorb  the 
heat  energy.  The  vapor  flows  to  the  rest  of  the  porous  surface  and 
releases  latent  heat  energy  as  it  condenses.  The  energy  is  dissipated 
into  space  or  to  the  environment  by  radiation  from  the  outer  surface. 
Because  of  the  high  latent  heat  of  the  working  fluid,  a  large  amount  of 
incoming  heat  can  be  absorbed  by  evaporation,  and  spread  to  the 
surrounding  surface  of  the  wall  to  be  dissipated  into  space,  without 
causing  the  temperature  in  the  wall  to  become  too  high  and  burn  out.  The 
positions  of  the  evaporator  and  condenser  are  not  fixed  nor  are  their 
areas.  These  depend  on  where  the  pipe  is  hit  by  the  incoming  heat  flux, 
and  how  large  the  surface  area  is  that  is  being  hit  by  the  heat  flux. 

Also,  this  kind  of  heat  pipe  has  no  adiabatic  section  (Fig.  3.1a).  The 

end  caps  have  been  removed  to  demonstrate  the  typical  interior  structure 
and  the  vapor  flow  pattern.  The  origin  of  the  coordinate  system  is  set 
at  the  center  of  the  heated  spot.  Figure  3.1b  shows  a  typical  wall 

temperature  profile  along  the  x  axis  at  y  =  0  and  the  vapor  temperature 

T  .  If  the  temperature  of  the  wall  is  higher  than  Tg,  it  serves  as  the 
evaporator;  if  the  wall  temperature  is  lower  than  T  ,  it  serves  as  the 

u 

condenser.  This  particular  phenomena  is  not  observed  for  conventional 
heat  pipes  operated  at  a  nearly  constant  wall  temperature. 

The  present  analysis  is  based  on  a  number  of  physical  assumptions. 


b  Typical  wa 


Ve  must  emphasize  the  wall  temperature  in  the  analysis  to  prevent  the 
heat  pipe  wall  from  melting.  The  heat  pipe  is  assumed  to  contain  only 
one  fluid  and  is  operating  under  steady- state  conditions.  This  means 
that  the  vapor  region  is  nearly  isothermal,  and  the  total  energy  input 
and  output  are  equal.  The  wick- liquid  matrix  is  assumed  to  be  very  thin, 
and  its  existence  is  represented  by  the  evaporation  heat  transfer 
coefficient,  Hg,  and  the  condensation  heat  transfer  coefficient,  H^. 
This  is  a  good  approach  for  the  purpose  of  designing  a  heat  pipe  because 
in  most  cases,  evaporation  and  condensation  heat  transfer  coefficients  in 
porous  media  can  only  be  obtained  from  experimental  data  (Dunn  and  Reay 
[1982];  Groll  et  al.  [1984];  Abhat  and  Seban  [1974]).  The  common 
measurable  parameters  in  experiments  on  heat  pipes  are  Tg,  Hg,  H^,  and 
the  wall  temperature.  The  results  of  the  analysis  based  on  the  above 
concept  can  be  readily  used  as  a  guide  when  designing  heat  pipes. 
Finally,  it  is  assumed  that  the  wall  thickness  is  thin,  and  the  radius  of 
curvature  of  the  surface  is  much  larger  than  h,  so  that  the  wall 
temperature  does  not  vary  substantially  with  radial  position  in  the  wall. 

In  the  Cartesian  coordinate  system  shown  in  Fig.  3.1a,  the 
two-dimensional  governing  equations  and  boundary  conditions  for  the  wall 
can  be  written  by  an  energy  balance  of  the  differential  elements  in  each 
section  of  the  heat  pipe.  For  the  evaporation  section,  except  the 
surface  beneath  the  heated  spot,  the  governing  equation  is 


d  r  d T-,  d  f  5T-, 

"o  k~o  +  “5  k  a 
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T  =  0 


(3.1) 
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where  T 


=  Tg  at  the  junction  of  the  evaporator  and  the  condenser 
segments,  e  is  the  emissivity,  S  is  the  thickness  of  the  wall,  Tq  is  the 
temperature  of  the  environment,  and  qg  is  the  evaporation  heat  flux.  Ve 
have  shown  that  the  evaporation  heat  transfer  coefficient  varies  with 
heat  flux  because  the  menisci  in  a  capillary  evaporator  recede  as  the 
heat  flux  increases.  An  examination  of  the  reported  data  shows  that  a 
power- law  boiling  relation  is  appropriate  for  relating  the  heat  flux  to 
the  evaporating  temperature  drop  in  a  heat  pipe  (Rosenfeld  [1987] ;  Dunn 
and  Reay  [1982]),  i.e.; 

%  =  »(T  -  Ts)b  (3.2) 

The  most  common  values  of  the  exponent  vary  from  1.0  to  2.0,  with 
that  of  liquid  metal  heat  pipes  remaining  1.0  and  the  coefficient  a  =  Hg 
in  many  cases  (Dunn  and  Reay  [1982];  Groll  et  al.  [1984];  Davis  and 
Ferrell  [1974]). 

For  the  wall  beneath  the  heated  spot,  the  corresponding  equation  is 
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where  q^  is  the  incoming  heat  flux  at  the  heated  spot. 

For  the  condenser  segment,  the  governing  equation  is 


d  r 

dx  kdx 


d  r  dT, 
dy  kdy 


(3.4) 


Vith  T  =  Tg  =  Tg  at  the  junction  of  the  evaporator  and  the  condenser 
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3T 

sections,  and  =  0  at  the  peripheral  boundary  of  the  heat  pipe  where  n 

is  the  normal  direction  of  the  periphery. 

The  vapor  temperature,  T  ,  is  an  unknown  in  equation  (3.4).  The 
fact  that  the  vapor  temperature  must  be  determined  requires  one  more 
equation  which  is  provided  by  an  energy  balance  over  the  entire  heat 
pipe.  A  steady- state  operation  has  no  energy  and  mass  accumulations,  and 
the  vapor  temperature  will  be  adjusted  according  to  the  heat  input  and 
the  ambient  conditions.  This  means  that  the  condensing  heat  transfer  is 
equal  to  that  of  evaporation,  and  all  the  heat  input  is  rejected  to  the 
ambient,  i . e . ; 

J JHe(T  -  Ts)di  =  JJHC(TS  -  T)dA  (3.5) 

AC 

Q„  =  |/e»(T4  -  T4)dA  (3.6) 

A 

where  Qjj  is  the  heat  input  through  the  spot;  Ag  and  A^  are  the  evaporator 
surface  area  and  the  condenser  surface  area  inside  the  heat  pipe, 
respectively;  A  is  the  outer  radiation  surface  area.  Note  that  A  is  not 
necessarily  equal  to  A^.  When  spot  heating  heat  pipes,  the  evaporator 
and  condenser  are  not  prescribed.  Usually,  the  heated  area  is  small,  and 
a  large  amount  of  heat  needs  to  be  spread  to  the  surrounding  surfaces  to 
be  dissipated.  Therefore,  A^  is  larger  than  the  spot  heating  area  and 
as  a  result,  A  >  A^,. 

3.3.2  Block  heated  heat  pipes  with  low  or  moderate  working  temperatures 
Block  heated  heat  pipes  are  normally  used  to  transport  energy  from 
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one  place  to  another  to  conserve  energy  or  to  cool  electronic  components. 
The  evaporator  and  condenser  segments  are  normally  separated,  which  is 
sir Jar  to  conventional  heat  pipes  except  for  the  localized  heating.  The 
working  temperature  is  comparatively  low  or  moderate,  which  means  that 
radiation  heat  transfer  is  not  important  in  the  analysis  of  these  heat 
pipes. 

Consider  the  heat  pipe  shown  in  Fig.  3.2.  It  has  an  evaporator 
section  of  length,  Lg,  a  wall  thickness,  6,  an  outside  radius,  R,  and  a 
block  heated  area  of  width,  Vjj ,  and  length,  L^. 

The  analysis  here  is  based  on  similar  assumptions  used  for  spot 
heated  heat  pipes  with  high  operating  temperatures.  With  a  thin  wall  and 
a  large  pipe  diameter,  the  problem  can  be  investigated  in  Cartesian 
coordinates. 

In  the  evaporator  section,  the  governing  equation  for  the  pipe  wall 
that  is  not  underneath  the  heated  block  is 


In  equation  (3.8),  we  have  employed  the  power- law  relation  for  the 
boiling  heat  flux.  In  this  situation,  b  is  greater  than  1.0.  The 
boundary  conditions  for  this  problem  are 
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(3.9) 


n 

^  =  0  at  x  =  0  and  x  =  Lg 


(3.10) 

Equation  (3.9)  defines  an  insulated  boundary  condition,  and  equation 
(3.10)  is  true  due  to  symmetry. 

Equations  (3.9)  and  (3.10)  can  be  nondimensionalized  with  the 
following  variables: 


(3.9a) 

I 

(3.10a) 

where  Cj  =  -a(xR)2Tgb  */5k,  =  (xR)2/Tg6k  and  Tg  is  the  reference 

temperature. 

The  above  model  is  justified  because  since  the  main  purpose  is  to 
transport  energy  from  one  place  to  another,  the  heat  dissipation  to  the  < 

ambient  from  the  evaporator  section  is  negligible.  Also,  because  of  the 
large  evaporation  heat  transfer  coefficient,  very  little  heat  is 
transported  from  the  evaporator  to  the  condenser  through  the  pipe  wall. 

Section  3.4  shows  the  numerical  results. 


T  =  (T  -  Ts)/Th,  X  =  x/xR,  Y  =  s/xR 
The  resulting  nondimensional  equations  are 


9  *  9  * 
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=  0 
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^  =  0  at  s  =  0  and  s  =  rR 
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3.4  Numerical  Results  and  Discussion 


The  governing  equations  for  spot  heated  and  block  heated  heat  pipes 
are  nonlinear  and  nonhomogeneous  and  require  an  iterative  procedure.  The 
problems  were  solved  by  using  the  f inite- difference  method  based  on  the 
control- volume  formulation  (Patankar  [1980]).  A  combination  of  the 
direct  method  (TDMA)  and  the  Gauss- Seidel  method  was  employed  to  solve 
the  discretization  equations.  In  some  special  cases,  under- relaxation 
parameters  were  used  to  control  the  advancement  of  the  solutions.  The 
grid  size  employed  in  the  program  varies  from  20  x  70  to  70  x  300 
depending  on  the  computational  domain.  To  determine  a  suitable  grid 
size,  the  computed  temperature  profiles  are  compared  for  a  number  of 
different  grid  sizes  for  the  same  problem.  The  maximum  errors  among 
these  grid  sizes  are  less  than  1.07,.  Also,  an  energy  balance  over  the 
entire  computational  domain  was  checked  for  every  computed  temperature 
field,  the  maximum  error  of  which  was  at  most  0.17.. 

Figure  3.3  shows  the  numerical  results  for  spot  heated  heat  pipes 
with  high  working  temperatures.  Like  conventional  heat  pipes,  the  heat 
input  has  a  strong  influence  on  the  heat  pipe  performance.  As  the  heat 
input  increases,  the  peak  wall  temperature  increases.  For  example,  the 
peak  temperature  with  =  8000  V  is  more  than  twice  as  high  as  that  with 
Qh  =  2000  V. 

In  conventional  heat  pipes,  the  area  available  for  heat  input  is 
comparatively  large,  the  input  heat  flux  is  comparatively  low,  and  the 
working  limitation  is  mainly  determined  by  the  total  energy  input.  Vith 
spot  heated  heat  pipes,  the  situation  is  different.  Figure  3.4  shows  the 
temperature  distribution  of  the  wall  with  a  very  high  incoming  heat  flux 
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(as  high  as  10  kV/cm  ).  Comparing  the  curve  of  H  =  Hg  =  =  10,000 
V/(m^-K)  and  Qg  =  2,750  V  in  Fig.  3.4  with  that  of  Qg  =  4000  V  in  Fig. 
3.3,  we  know  that  even  though  the  heat  input  Qg  =  2750  V  in  Fig.  3.4  is 
less  than  that  of  Qg  =  4000  V  in  Fig.  3.3,  its  peak  temperature  is  still 
several  hundred  degrees  higher  because  of  the  smaller  area  that  is 
heated.  On  the  other  hand,  Fig.  3.4  shows  that  a  higher  heat  transfer 
coefficient  can  reduce  the  peak  temperature,  as  is  expected. 

Also  shown  in  Fig.  3.4  is  the  influence  of  Ag  on  the  peak  wall 
temperature.  The  curve  with  a  much  lower  peak  wall  temperature  in  Fig. 
3.4  is  subject  to  the  same  heat  flux  as  that  of  the  other  three  curves 
(q^,  =  10  kV/cm  ),  but  a  smaller  heating  area  (Qg  is  also  smaller).  The 
resultant  peak  wall  temperature  is  almost  1,000  degrees  lower  than  that 
with  the  larger  heating  area.  In  this  study,  the  heating  area  is  square. 
Because  of  the  small  heating  areas,  its  shape  is  less  important. 

Unlike  conventional  heat  pipes,  the  thermal  conductivity  of  the  wall 
influences  the  wall  temperature  distribution  greatly  for  spot  heated  heat 
pipes,  especially  in  the  case  of  a  very  high  incoming  heat  flux,  as  shown 
in  Fig.  3.5.  When  the  thermal  conductivity  of  the  wall  is  small,  the 
peak  temperature  would  jump  intolerably  high  and  the  surface  would  be 
burned  out.  On  the  other  hand,  with  a  large  wall  thermal  conductivity, 
the  peak  temperature  decreases  sharply.  This  is  not  surprising  because  a 
large  amount  of  heat  needs  to  be  transferred  to  the  surrounding  wall 
through  a  very  small  area  by  conduction. 

In  conventional  heat  pipes,  better  cooling  conditions  and  a  larger 
cooling  surface  can  ameliorate  the  performance  and  result  in  a  lower 
working  temperature.  This  is  also  true  for  spot  heated  heat  pipes.  The 
trend  is  well  illustrated  in  Fig.  3.6,  where  a  larger  surface  area 
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(radiation  heat  transfer  area)  corresponds  to  a  lower  working  temperature 
T  .  Also,  the  peak  wall  temperature  is  decreased  accordingly.  But,  if 
we  pay  attention  to  the  area  around  the  center  of  the  heated  spot,  we  may 
notice  that  the  temperature  difference  between  the  peak  temperature  and 
the  working  temperature  for  different  curves  almost  remains  the  same. 
This  local  phenomenon  obviously  results  from  the  localized  heat  input 
characteristics  of  the  heat  pipe.  The  most  important  factor  which 
determines  the  performance  is  the  working  conditions  at  the  heated  spot. 

To  estimate  the  validity  of  spot  heating  heat  pipes  to  reduce  the 
wall  temperature,  we  compared  the  wall  temperature  of  a  plate  that  is  not 
a  heat  pipe  with  that  of  a  plate  heat  pipe.  The  curve  with  the  solid 
circle  legend  in  Fig.  3.7  is  the  wall  temperature  profile  of  a  spot 
heated  heat  pipe  under  normal  working  conditions,  while  the  curve 
indicated  with  H  =  0  is  the  temperature  profile  of  a  plate  that  is  not  a 
heat  pipe  (H^  =  =  H  =  0),  with  other  conditions  being  the  same. 
Obviously,  the  peak  wall  temperature  of  the  surface  adopting  heat  pipe 
technology  is  reduced  significantly. 

Vhen  we  compare  the  conventional  heat  pipes,  it  is  clear  that  spot 
heated  heat  pipes  have  small  evaporator  surfaces,  very  high  evaporation 
heat  fluxes,  large  condenser  surfaces  and  vapor  passages.  Because  of 
these  factors,  the  boiling  limit  becomes  the  most  important  operating 
limit  of  these  heat  pipes.  In  Fig.  3.7,  the  curve  indicated  with 
ATc  =  100  K  assumes  that  the  boiling  limitation  is  reached  for  the  local 
evaporating  surface  when  T  -  T  >  AT  =  100  K.  Vhen  this  occurs,  the 
porous  media  at  that  point  is  assumed  to  be  completely  dry,  and  no 
evaporation  takes  place.  For  the  curve  with  the  solid  circle  legend, 
this  restriction  has  not  been  imposed  on  AT,  and  the  pipe  works  properly 
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,3.7  TEMPERATURE  PRORLES  WITH  DIFFERENT  WORKING  CONDITIONS 


with  the  other  conditions  being  the  same.  Once  the  boiling  limitation  is 
reached  and  evaporating  surface  becomes  dry,  the  wall  temperature  at  the 
heated  spot  will  jump  thousands  of  degrees  higher  than  that  of  heat 
pipes  under  normal  working  conditions.  From  the  viewpoint  of  heat  pipe 
design,  most  difficulties  arise  from  the  avoidance  of  the  incipience  of 
boiling  in  the  porous  wick  of  the  pipe.  Special  care  must  be  taken  to 
properly  design  the  structure  of  the  porous  media,  to  choose  suitable 
working  fluids,  and  to  insure  the  wettability  of  the  wick  and  the  wall  to 
increase  the  boiling  limitation. 

In  addition,  the  emissivity  of  the  surface  and  the  thickness  of  the 
wall  have  strong  effects  on  the  temperature  profile  of  the  wall. 
Increasing  the  emissivity  will  reduce  the  working  temperature,  and 
therefore  the  peak  wall  temperature.  The  value  of  e  in  this  numerical 
calculation  was  taken  as  0.8.  Also,  a  larger  thickness  of  the  wall  will 
alleviate  the  temperature  jumps  at  the  center  of  the  heated  spot,  but 
this  is  not  practical  in  many  applications. 

Figures  3.8  through  3.10  show  the  numerical  results  for  heat  pipes 
with  localized  heat  input  working  under  low  or  moderate  temperatures.  In 
this  situation,  the  heated  area  is  comparatively  large,  and  the  heat  flux 
is  comparatively  low,  so  that  the  temperature  jump  is  not  so  severe  as 
that  for  spot  heating  heat  pipes  with  high  working  temperatures. 

Figure  3.8  shows  the  comparison  of  the  numerical  results  of  the 
circumferential  wall  temperature  profile  with  the  experimental  data  from 
the  paper  by  Rosenfeld  [1987].  The  experiment  was  conducted  with  a 
narrow  line  heater  at  the  evaporation  section  of  the  heat  pipe.  The 
evaporation  heat  flux  relation  is  also  taken  from  that  experiment,  with 


137 


Qh=950  (W) 
(5=0.003175  (m) 
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«E  =  a(T  -  Ts)b 

where  a  =  17330,  and  b  =  1.215. 

The  agreement  between  the  numerical  results  and  the  experimental 
data  is  excellent.  The  line  width  surely  influences  the  performance  of 
the  heat  pipe.  With  the  width  of  the  heated  line  becoming  larger,  the 
temperature  distribution  along  the  circumferential  direction  becomes 
smoother,  as  shown  in  Fig.  3.9.  Among  the  different  heating  widths,  half 
heating  (sg  =  Vg/rR  =  0.5)  is  of  special  interest  in  many  applications. 
Obviously  with  a  uniform  input  heat  flux  and  a  large  evaporation  heat 
transfer  coefficient,  the  temperature  profile  of  the  wall  beneath  the 
heated  block  is  nearly  smooth,  and  the  working  conditions  of  this  half  of 
the  evaporator  are  nearly  the  same  as  that  of  heat  pipes  with  a  uniform 
heat  input. 

Figure  3.10  shows  the  performance  of  block  heated  heat  pipes  for 
different  values  of  the  wall  thermal  conductivity.  Unlike  the  spot 
heated  heat  pipes  shown  in  Fig.  3.5,  the  thermal  conductivity  of  the  wall 
has  little  effect  on  the  temperature  distribution.  The  reason  is  that, 
because  of  the  large  boiling  heat  transfer  coefficient,  most  of  the  input 
heat  was  absorbed  by  the  evaporating  surface  beneath  the  heated  block, 
and  only  a  small  amount  of  heat  is  spread  to  the  wall  that  is  not  heated. 
This  is  more  pronounced  as  the  wall  thermal  conductivity  decreases.  For 
the  working  conditions  indicated  in  Fig.  3.10,  with  k  =  394  V/(m-K), 
86.27,  of  the  heat  passes  through  the  wall  under  the  heated  block,  while 
with  k  =  45  V/(m-K),  967.  of  the  heat  passes  through  the  wall  under  the 
heated  block. 
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3.5  Conclusions  and  Remarks 


1.  The  use  of  a  heat  pipe  is  an  excellent  method  to  protect  a  surface 
from  burning  out  when  the  surface  is  spot  heated.  The  peak  wall 
temperature  is  greatly  reduced  due  to  the  operation  of  the  heat  pipe. 
The  parameters  which  influence  the  performance  are  Qg,  k,  Hg,  H^,,  e 
and  <5.  For  a  fixed  heat  input  larger  values  of  k,  A^,  Hg,  H^,  or  6 
can  reduce  the  peak  wall  temperature,  and  larger  values  of  e  and  the  heat 
pipe  surface  area  result  in  a  lower  working  temperature.  Because  of  the 
localized  heating  characteristics  of  heat  pipes,  a  temperature  jump  is 
inevitable  at  the  center  of  the  heated  spot  and  results  in  a  large  AT  and 
a  high  at  the  heated  location.  Special  measures  must  be  taken  to 
increase  the  boiling  limitation  of  the  heat  pipe.  Otherwise,  the 
temperature  at  the  heated  location  will  jump  intolerably  high.  More  work 
needs  to  be  done  on  the  structure  of  the  porous  media,  the  wettability  of 
the  wick  and  wall,  and  the  vapor  flow  pattern  in  the  pipe  for  this 
special  kind  of  heat  pipe. 

2.  The  numerical  results  for  localized  heat  input  heat  pipes  working 
under  low  or  moderate  temperatures  agree  well  with  the  existing 
experimental  data  and  can  be  used  to  improve  the  prediction  of  heat  pipe 
performance  under  localized  heating.  With  a  large  evaporation  heat 
transfer  coefficient,  most  of  the  input  heat  passes  through  the  wall 
under  the  heated  block. 


142 


IV.  A  NUMERICAL  ANALYSIS  OF  STEFAN  PROBLEMS  FCR 
GENERALIZED  MULTI- DIMENSIONAL  PHASE- CHANGE  STRUCTURES 
USING  THE  ENTHALPY  TRANSFORMING  MODEL 


4.1  Summary 


An  enthalpy  transforming  scheme  is  proposed  to  convert  the  energy 
equation  into  a  nonlinear  equation  with  the  enthalpy,  E,  being  the  single 
dependent  variable.  The  existing  control- volume  finite  difference 
approach  has  been  modified  to  apply  it  to  the  numerical  performance  of 
Stefan  problems.  The  model  has  been  tested  by  applying  it  to  a 
three-dimensional  freezing  problem.  The  numerical  results  are  in  good 
agreement  with  those  existing  in  the  literature.  The  model  and  its 
algorithm  are  further  applied  to  a  three-dimensional  moving  heat  source 
problem  showing  that  the  methodology  is  capable  of  handling  complicated 
phase- change  problems  with  fixed  grids. 


4.2  Introduction 


Heat  flow  and  diffusion  with  melting  and  solidification  are  of  great 
importance  in  many  industrial  applications.  Examples  are  casting, 
welding,  thermal  energy  storage  units,  heat  pipe  start-up  from  frozen 
state,  etc.  The  last  two  operations  were  the  major  motivation  for  the 
present  study.  Phase- change  processes  may  produce  solid  and  liquid  phase 
regions  which  have  extremely  complex  appearances.  Also,  it  is  not 
possible  to  predict  a  priori  what  the  phase- change  front  evolving  in  time 
will  look  like.  Therefore,  exact  analytical  solutions  for  these  types  of 
nonlinear  problems  are  available  only  for  some  simplified  and  idealized 
systems.  Apparently,  numerical  methods  are  the  only  practical  method  to 
handle  the  general  melting  and  freezing  problems  providing  that  we  can 
successfully  trace  the  moving  interface. 

The  numerical  methods  used  to  solve  phase- change  problems  might  be 
divided  into  two  main  groups.  The  first  group  is  called  strong  numerical 
solutions.  The  focus  here  is  on  applying  f inite- difference  techniques  to 
the  strong  formulation  of  the  pro  ess,  locating  fronts  and  finding 
temperature  distributions  at  each  time  step  or  employing  a  transformed 
coordinate  system  to  immobilize  the  moving  interfaces  (Okada  [1984],  Ho 
and  Chen  [1986]).  These  methods  are  applicable  to  those  processes 
involving  one  or  two  phases  in  one  space  dimension  which,  with  the  use  of 
cumbersome  schemes,  are  being  applied  to  two-dimensional  cases  as  well. 

The  second  group  is  called  weak  numerical  solutions  (Shamsundar  and 
Sparrow  [1975],  Crowley  [1978],  Voller  and  Cross  [1981],  Hsiao  and  Chung 


144 


[1984],  Keung  [1980],  and  Hsiao  [1984]),  These  methods  allow  ns  to  avoid 
paying  explicit  attention  to  the  nature  of  the  phase  change  front.  They 
appear  to  have  great  flexibility  and  are  easily  extended  to 
multi- dimensional  problems.  In  this  group,  the  most  important  and  widely 
used  method  is  the  enthalpy  method.  The  advantages  of  the  enthalpy 
reformulation  are  that  the  problem  to  be  solved  is  formulated  in  a  fixed 
region,  and  no  modification  of  the  numerical  scheme  is  necessary  to 
satisfy  the  conditions  at  the  moving  phase- change  interface. 
Furthermore,  this  method  is  especially  suitable  both  for  the  problems 
where  the  phase  change  occurs  at  a  single  temperature  and  the  problems 
where  the  phase  change  occurs  over  a  temperature  range. 

Most  of  the  previous  enthalpy  models  usually  treated  the  enthalpy  as 
a  dependent  variable  in  addition  to  the  temperature  and  discretized  the 
energy  equation  into  a  set  of  equations  which  contain  both  E  and  T.  For 
the  implicit  schemes,  they  actually  treated  all  of  the  terms  containing  T 
-  T(E)  as  a  constant  heat  source  term  in  the  energy  equation  during 
iterations  at  each  time  step.  This  may  cause  some  problems  for 
convergence  when  T  =  T(E)  is  complicated  and  physical  properties  change 
significantly  as  is  the  case  of  frozen  heat  pipe  start-up,  or  when 
boundary  conditions  are  severe.  Furthermore,  when  the  energy  equation 
contains  a  convective  term,  the  previous  methods  have  difficulties  in 
handling  the  relationship  between  the  convective  and  diffusive  terms 
because  of  the  two  dependent  variable  nature  of  the  equation. 

In  this  section,  a  simple  strategy  is  proposed  to  transform  the 
energy  equation  into  a  nonlinear  equation  with  a  single  dependent 
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variable  E.  Thus,  solving  a  phase- change  problem  is  equivalent  to 
solving  a  nonlinear  enthalpy  equation,  and  existing  algorithms  are 
readily  applicable  with  some  modifications. 
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4.3  Enthalpy  Transformation  of  the  Energy  Equation 


The  energy  equation  governing  three-dimensional  laminar  flow  with  no 
viscous  dissipation  and  with  incorporating  the  continuity  equation  in  the 
Cartesian  coordinate  system  (Kays  and  Crawford  [1980],  Patankar  [1980]) 
is 


3(pE)  d  d  d  d  91  d  91  d  91 

~Si  *  Sx  (',uE)  *  Sj  OvE>  *  Si  ^  =  Si  (kSi]  *  J^kSy)  +  Si  (kSi> 

(4.1) 


with  the  state  equation 


dE 

=  C(T)  (4-2) 

In  the  case  of  constant  specific  heats  for  each  phase,  and  that  the  phase 
change  occurs  at  a  single  temperature,  we  have 


[  T  +  E/c  E  <  0 

I  m  '  s 

T  0  <  E  <  H 

m 


L  Tm  -  (E  -  H)/C/  E  > 


H 


(solid  phase) 
(mushy  phase) 
(liquid  phase) 


(4.3) 


where  is  the  melting  temperature.  In  the  above  relation,  we  have 
selected  E  =  0  to  correspond  to  phase  change  materials  in  their  solid 
state  at  temperature  T  . 
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The  "Kirchoff"  temperature  (Solomon  et  al.  [1966])  is  introduced  as 


follows 


MT  -  V 


T+  =  r  k(if)d» I  =  0 

T_ 


T  <  T 


T  =  T 


(4.4) 


k^(T  -  Tm)  T  >  T„ 


Transforming  eq.  (4.3)  with  the  definition  given  in  eq.  (4.4)  results  in 


ksE/cs 


E  <  0 


k£(E  -  H)/C/ 


0  <  E  <  H 


E  >  H 


(4.5) 


and  eq.  (4.1)  becomes 


d(pE)  +  d 


d2j*  d*T* 


di  (PUE)  ’  ^  (pvE)  T  (?vE)  "  ^  *  ^2  '  3J5 


(4.6) 


Now,  let  us  introduce  an  enthalpy  function  as  follows 


T  =  r(E)E  +  S(E) 


(4.7) 


For  the  phase  change  occurring  at  a  single  temperature,  we  have 


k  /  c 
s'  s 


E  <  0 


T(E)  =  0 


0  <  E  <  H 


(4.8) 


E  >  H 
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and 


r  o 


E  <  0 


S(E) 


0 


0  <  E  <  H 


E  >  H 


(4.9) 


Upon  substituting  eq.  (4.7)  into  eq.  (4.6)  and  noticing  that,  for 
example, 


dn  d  p(rE+s)' 

dx 2  dx  dx 


d 2  32S 

=  ^2(rE)  + 


we  have 


fl(/>E)  d  d  d  =  52TE  +  02rE  +  02  TE  +  p 

”3T~  +  3*  ((,uE)  *  ^  (pvE)  *  S  (/,wE>  "sr^Tp 


(4.10) 


with 


P  = 


<92S 

0X2 


02  s 

0y2 


02S 

0Z2 


r  -  r(E)  , 


and  S  =  S(E) . 


The  energy  equation  has  been  transformed  into  a  nonlinear  equation 
with  single  dependent  E.  The  nonlinearity  of  the  phase- change  problem  is 
evident  in  the  above  equation. 


In  the  liquid  region  away  from  the  moving  front  as  indicated  in  the 
numerical  domain  of  Fig.  4.1,  eq.  (4.10)  reduces  to  the  normal  linear 
energy  equation  as 
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d  d  d  d  91  d  dl  d 

“3T“  +  +  Ty^ly^  +  ^WwE)  "  IS)  +  Ty^i  Ty">  +  Tz^l 

ffl 

dz )  (4-11) 

Also,  in  the  solid  region  eq.  (4.10)  reduces  to 


^(PSE) 

dt 


+  J^gUE) 
dx  s 


+  _f(/>gVE) 

dy  s 


d(  v v  <?,,  dTx 
+  _(^qwE)  =  _(k  _  ) 


dz 


dx  dx 


d(.  dT,  d(,  <?Tx 
dy  sdy  dz  sdz 
(4.12) 


In  the  moving  front  region  (the  region  between  the  two  dashed  lines 
as  indicated  in  Fig.  (4.1),  eq.  (4.10)  is  nonlinear.  This  agrees  with 
the  well-known  fact  that  the  nonlinearity  of  phase- change  problems  is  due 
to  the  existence  of  a  moving  interface  (Ozisik  [1980]). 

The  method  proposed  is  not  restricted  to  the  forms  for  T(E)  and  S(E) 
given  by  eq.  (4.8)  and  eq.  (4.9).  With  different  conditions  and 
assumptions,  they  have  different  expressions.  For  example,  if  phase 
changes  occur  over  a  temperature  range  (such  as  alloys),  as  shown  in  Fig. 
4.2,  with  constant  specific  heats  for  each  phase,  we  have 

E/c  E<0  (solid  phase) 

Q 

ATE/ ( H  +  cmAT)  0<E<H+cmAT  (mushy  phase)  (4.13) 
E/c^-  [H+(cm- c^)AT] /c^  E>H+cmAT  (liquid  phase) 

Here,  we  have  selected  E  =  0  to  correspond  to  phase- change  materials 
in  their  solid  state  at  temperature  T^.  Then  T^  =  is  defined 
as  the  melting  temperature,  AT  =  T^  is  the  melting  temperature 
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range,  and  Cm  is  the  specific  heat  for  the  mushy  phase. 


The  "Kirchoff"  temperature  is  introduced  as 


T+  =  (1  k(i |)  dij  = 
*1 


'  ks  (T-Tj) 

K„  (I-tp 

(T-T,) 


T  <  T, 

T,  <  T  <  T2 


(4.14) 


where  km  is  the  thermal  conductivity  for  the  mushy  phase.  The 
transformation  procedure  is  the  same  as  that  of  phase  change  at  a  single 
temperature,  and  the  resulting  equation  is  still  eq.  (4.10)  with 
different  expressions  for  T(E)  and  S(E). 
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(4.15) 


(4.16) 


In  the  above  relations  for  the  mushy  region,  a  linear  change  was 
assumed.  In  real  systems,  they  may  take  more  complicated  forms. 
However,  this  is  outside  the  scope  of  this  report. 
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4.4  Numerical  Scheie 


4.4.1  Phase  change  without  convective  terms 

To  demonstrate  the  methodology,  let  us  consider  a  phase- change 
problem  in  one  space  dimension.  In  this  case,  eq.  (4.10)  reduces  to 


5t 

:  92  (TE) 
ox2 

+  52s 

5x2 

(4.17) 

with 

r  =  r(E) 

and  S  =  S(E) 

The  discretization  of  the  above  equation  employs  the  control- volume 
finite- difference  approach  described  by  Patankar  [1980].  In  this 
methodology,  the  discretization  equations  are  obtained  by  applying 
conservation  laws  over  finite  size  control  volumes  surrounding  the  grid 
nodes,  and  integrating  the  equation  over  the  control  volumes,  i.e. 

///  p  9E  dV  =  ///  (  52(rE)  +  52S  )  dV  (4.18) 
AV  5t  AV  5x2  5x2 


Using  a  fully  implicit  scheme  and  referring  to  Fig.  4.3,  we  have 


J  IjP—  dV  =  P^x 

AV  5t 


Ep  - 
At 


(4.19) 
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5rE  x  ,  5rE  X  rEEE  '  rPEP  rPEP  '  rvEv 

"35  6  "E  tf  " 


(4.20) 


|ff 

jli  5x2  ^ 


(4.21) 


154 


w 

V 

e 

VI 

// 

// 

45#' 

e: 

-  AX  - -j 

1 

FIG.  4.3  GRID- POINT  CLUSTER  FOR  THE 
ONE- DIMENSIONAL  PROBLEM 


155 


Thus 


apEp  -  a-gEg  +  &yEy  +  b  (4.22) 

o 

with  Ep  denoting  the  old  value  of  E  at  grid  point  P 
rE  FV 

“  _  5  —  _ 

Ux)e  Ux)w 

and 

a  _  FP  rP  +  pAx 

P  W~e  "At 

4.4.2  Phase  change  with  convective  terms 

In  this  case,  a  one  space  dimensional  problem  will  also  be 
considered  as  a  demonstrative  example.  The  governing  equation  is 

dpi  d  _  d*  52S  (4.23) 

~dt  cbc  W  (EE)  W 

Since  the  total  flux  in  the  above  equation  J  =  pul  -  ^(TE)  is  different 
from  the  conventional  total  flux  J  =  pu§  -  the  usual  method  to 

obtain  the  convection- diffusion  expression  is  not  applicable.  Also,  the 
coefficient  T  is  small  in  most  cases.  To  handle  convection- diffusion 
situations  and  ensure  physically  realistic  solutions,  we  employed  a 
scheme  similar  to  the  upwind  scheme.  The  discretization  equation  is 
written  as 

apEp  -  a^Eg  +  ayEy  +  b  (4.24) 


pAx ED°  Sr  -  SD  SD  -  Sir 

b  =  _ +  E  P  -  P  V  ,  (4.22a) 

At  (<5x)e  (Ax)w 

(4.22b) 
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where 


_  o 

ap  '  aPE  +  aPV  +  At  p P 

aE  =  I'EDe  +  max  [-Fe,  0],  apE  =  rpDe  +  max  [-Fe,  0] 

ay  =  FyD^  +  max  [F^,  0] ,  apy  -  FpD^  +  max  i_F ^ ,  0] 

AxE£  SE  -  Sp  Sp  -  Sy 
St~  +  (6x)^  U*)w 

De  =  1/(^x)e’  Fe  =  (He>  Dw  =  1/(Sx)v>  fw  =  ^u)w 

The  greater  of  a  and  b  is  given  by  max  [a,  b] ,  Ep  denotes  the  old 
value  of  E  at  grid  point  P,  and  />p  denotes  the  old  value  of  p  at  grid 
point  P.  Clearly,  no  special  treatment  is  needed  for  solving  the 
velocity  using  momentum  equations. 


4.4.3  Phase  change  for  multi- dimensional  problems 

Having  described  the  discretization  equation  for  one  space 
dimensional  problems,  we  can  now  write  a  discretization  equation  based  on 
the  general  differential  equation  (4.10)  for  multi  dimensional  problem.-,, 
with  E,  V,  N,  S,  T,  and  B  representing  the  "east,"  "west,"  "north," 
"south,"  "top,"  and  "bottom"  neighbors  of  node  P,  respectively.  The 
corresponding  discretization  equation  is 

apEp  =  aEEE  *  ayEy  +  a^Ej,  +  agEs  +  aTET  *  apEB  +  b  (4.25) 


w  h  e  r  e 

Ax  Ay  Az  o 

aP  aPE  *  aPV  '  aPN  ‘  aPS  ‘  aPB  +  aPT  f  At  --p? 
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aE  "  Ve  +  max  f'Fe’  °] »  aPE  =  FPDe  +  max  f'Fe’  °] 

aV  =  FVDw  +  max  fFw’  °1’  aPV  =  FPDw  +  max  tFw’  °] 

aN  =  FNDn  +  max  t'Fn’  °3  ’  aPN  =  FPDn  +  max  t'Fn’  °] 

aS  =  PSDS  +  max  [Fg,  0],  apg  =  rpDg  +  max  [Fg,  0] 

aT  =  FT^t  +  max  ^Ft’  aPT  =  FP^t  +  max  ^Ft’  ^ 

aB  =  rBDb  +  max  [Fb,  0],  apfi  =  rpDb  +  max  [Fb,  0] 

AxAyAzpp 

b  =  — ^ -  EP  +  MW  "  MW  +  WV 

DS(VSS)  +  ®t^T  ^P^  bb^P~EB) 

The  flow  rates  and  conductances  are  defined  as 


Fe  =  (/»0e  AyAz 

Fw  =  (Hh  AyAz 
Fn  =  (Hn  AzAx 

Fs  =  ^s  AzAx 
Ft  =  (Ht  AxAy 

Fb  =  ^b  AxAy 


D  =  Ayaz 

e  ^ 
D  _  AyAz 
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p  _  AzAx 

'mw, 

p  _  AzAx 

5  =  wj; 

D,  = 

1  to; 

i),  s*4 y 

b  TOt 


(4.26) 


Because  of  the  nonlinearity  of  the  above  equation  and  the  implicit 
nature  of  the  scheme,  iterations  are  needed  at  each  time  step.  This 
procedure  is  the  same  as  that  which  solves  a  nonlinear  equation,  and  is 
as  follows: 


158 


r 

1)  Let  E  represent  the  E  field  as  it  exists  at  the  beginning  of  the 
kth  iteration. 

2)  From  these  values,  calculate  tentative  values  of  F  and  S  according 
to  their  relations  with  E,  using  eqs.  (4.8)  and  (4.9),  or  eqs. 
(4.15)  and  (4.16). 

3)  Solve  the  nominally  linear  set  of  discretization  equations  to  get 
new  values  of  E^+^. 

4)  Return  to  step  1  and  repeat  the  process  until  further  iterations 
cease  to  produce  any  significant  change  in  the  values  of  E. 


4.5  Application  of  the  Methodology  to  the  Example  Problems 


To  demonstrate  the  present  scheme,  the  proposed  methodology  has  been 
applied  to  two  separate  phase- change  problems.  The  first  is  a 
three-dimensional  freezing  problem,  and  the  other  is  a  three-dimensional 
moving  heat  source  problem.  In  the  two  problems  considered  herein,  the 
thermal  physical  properties  such  as  k  and  C  are  assumed  to  be  constant  in 
each  phase  but  may  differ  among  the  solid,  mushy  and  liquid  phases,  while 
the  density,  p,  is  considered  the  same  for  each  phase. 

4.5.1  The  three-dimensional  freezing  problem 

Consider  a  liquid  initially  at  its  melting  temperature,  T  ,  in  a  bar 
with  a  uniform  square  cross  section  and  adiabatic  ends  as  shown  in  Fig. 
4.4a.  The  surface  is  suddenly  exposed  to  a  uniform  wall  temperature 
below  the  fusion  temperature  and  freezing  takes  place  immediately. 
Because  of  the  symmetry  of  the  geometry,  only  a  quarter  of  the  bar  is 
considered  as  shown  in  Fig.  4.4b.  To  facilitate  comparison,  the 
dimensionless  parameters  are  chosen  to  be  the  same  as  those  used  by  Hsiao 
and  Chung  [1984] ,  i ,e. 

0,  =  (T.  -  T  )  /  (T  -  T  )  =  1  and  St  =  c  (T  -  T  )  /  H  =  1. 

At  the  middle  plane  of  the  bar  in  the  z-direction,  the  temperature 
distribution  is  two-dimensional.  Figure  4.5  shows  the  interface  position 
as  a  function  of  time  along  the  diagonal  for  the  present 
three-dimensional  modeling.  The  two-dimensional  results  given  by  Hsiao 
and  Chung  [1984]  using  the  equivalent  heat  capacity  model  and  given  by 
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(a)  Bar  of  liquid  with  an  uniform  square  cross  section 


O  T 
1  w 


KH 

(b)  One- quarter  of  the  bar  used  for  the  computational 
domain  due  to  the  symmetry  of  the  problem. 
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Crowley  [1978]  using  the  enthalpy  model  are  also  included  in  the  same 
figure.  The  agreement  among  these  solutions  is  excellent. 

Consider  the  same  problem  with  different  initial  conditions  and 
physical  properties.  The  dimensionless  parameters  are  k ^/k g  =  0.9, 
at/as  --  0.9,  0,  --  (T;  -  Ty)/(Tm  -  Ty)  =  9/7,  and  St  ,  cs(Tk  Ty)/ll  =  2, 
Figure  4.6  shows  the  interface  position  as  a  function  of  time  along  the 
diagonal.  Also  included  in  Fig.  4.6  are  solutions  obtained  by  Hsiao  and 
Chung  [1984],  and  by  Keung  [1980].  Again,  the  present  three-dimensional 
solution  agrees  well  with  the  results  of  those  two-dimensional  studies. 

Ve  base  the  above  calculation  on  the  phase  change  which  occurs  at  a 
single  temperature.  Assuming  that  the  phase  change  takes  place  over  a 
temperature  range  of  AT  =  20  K,  the  same  calculation  is  conducted,  as 
also  shown  in  Fig.  4.6.  Plainly,  the  present  model  is  insensitive  to  the 
phase- change  temperature  range.  If  the  temperature  range  is  small 
enough,  we  can  expect  the  same  result  as  that  of  the  single  temperature 
case.  This  is  the  case  for  the  present  model.  The  calculation  is 
conducted  with  AT  -  2  K,  and  the  result  is  almost  identical  with  that  of 
the  single  temperature  case  in  Fig.  4.6. 

The  grid  size  employed  in  the  above  two  cases  is  20  x  20  x  30,  and 
the  discretization  equations  are  solved  by  the  Gauss- Seidel  method.  The 
physical  properties  of  the  mushy  phase  are  taken  as  the  average  of  those 
of  the  solid  and  liquid  phases.  The  time  step  limit  is  not  encountered 
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in  the  calculations.  The  dimensionless  time  steps  can  be  on  the  order  of 
0.1,  and  real  time  steps  can  be  as  large  as  several  days. 


4.5.2  Three-dimensional  phase- change  problem  with  moving  heat  source 
As  indicated  in  Fig.  4.7,  a  source  of  heat  moves  over  the  surface  of 
the  plate  with  speed  U.  Because  of  intense  heating,  the  material  under 
the  heat  source  melts.  It  is  important  to  determine  the  molten  depth  for 
the  given  velocity,  heat  source  power  and  its  diameter,  as  well  as  the 
material  properties. 


With  the  coordinates  fixed  at  the  center  of  the  moving  heat  source, 
eq.  (4.10)  is  applicable  for  this  problem.  To  simulate  the  circular  heat 
source,  the  equation  has  been  transformed  into  the  form  for  the 
cylindrical- polar  coordinate  system. 

The  governing  equation  is 


«(E/>)  +  1  S<rVE^  1  a<v/E) 
dt  v  dr  r  dO 


d  ra(rE)' 
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p  =  1  d 
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vf  =  -  Ucos#  ,  =  UsinO 


165 


MOLTEN 
PHASE - 


U 

Q 


x 


(a)  Side-View  of  the  Computational  Domain  with  the 
Heat  Source  Moving  at  Speed  U. 


(b)  Top- View  of  the  Computational  Domain  with  the  Heat  source 
stationary  and  the  Computational  Domain  moving  at  speed  -U. 
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r  =  T(E)  ,  S  =  S(E)  are  given  by  eqs.  (4.8)  and  (4.9), 
respectively. 

The  initial  and  boundary  conditions  are  as  follows 


E.  =  c  (T.  -  T  ) 

1  S  v  1 

t  =  0 

q  =  qh 

z  =  0  ,  r  <  R^ 

^2 

II 

O 

z  =  0  ,  r  >  R^ 

li 

o 

z  -6 

R  -  Ei 

r  =  R0  ,  -  90°  <  0  <  90° 

<5(rE  ♦  S)  =  0 
dr 

r  =  R0  ,  90°  <  0  <  270° 

The  radius  R0  must  be  sufficiently  large  such  that  the  region  r  > 
o  o 

R0 ,  -90  <  0  <  90  is  unaffected  by  the  moving  heat  source.  Also,  the 

last  boundary  condition  implies  that  the  upwind  scheme  is  used,  and  the 
diffusive  term  is  neglected  for  the  outflow  boundary.  The  calculation 
proceeds  with  grid  size  32  x  50  x  12  and  time  step  0.1  second.  Other 
parameters  are 

U  =  0.3  m/s,  ajap  =  1.44,  St  =  c.(T  T - ) / II  =  0.126 
s  c  s  m  l 

R^  =  0.005  m,  R0  -  0.25  m,  (j  -  11.80  kV 


167 


Figure  4.8  shows  isotherms  of  the  dimensionless  temperature  0^  =  (T 
-  T._) /  (Tm  -  T^)  for  t  =  0.1  second  at  the  x  -  x  plane  indicated  in  Fig. 
4.7  (b)  (i.e.,  9  =  0°,  and  9  =  180°).  The  center  of  the  heat  source  is 
located  at  x  =  0.0.  The  solid  line  labelled  with  0  indicates  the  melting 
front  at  this  time,  while  the  dashed  line  labelled  with  -1  is  a  boundary 
beyond  which  the  temperature  field  is  unaffected  by  the  moving  heat 
souice.  After  about  0.5  second,  the  steady- state  condition  is  reached. 
Figure  4.9  shows  the  steady- state  isotherms  of  the  dimensionless 
temperature  at  the  same  plane.  The  melting  front  line  becomes  flat  in 
the  portion  of  x  <  0. 
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FIG.4.8  ISOTHERMS  OF  THE  SOLUTION  FOR  t=0.1s  AT  X-X  PLANE 


L9  STEADY  STATE  ISOTHERMS  OF  THE  SOLUTION  AT  X-X  PLANE 


4.6  Conclusions 


The  enthalpy  transforming  model  proposed  in  this  section  proves  to 
be  capable  of  handling  complicated  phase- change  problems  occurring  both 
at  a  single  temperature  and  a  temperature  range  with  fixed  grids.  Due  to 
the  one  dependent  variable  nature  of  the  transformed  equation,  the 
convection  and  diffusion  situations  can  be  handled  with  appropriate 
algorithms.  Ve  have  compared  the  numerical  results  existing  in  the 
literature  with  good  agreement,  showing  that  the  present  model  can 
properly  predict  the  phase- change  processes.  The  advantage  of  this  model 
based  on  enthalpy  is  that  it  allows  us  to  avoid  paying  explicit  attention 
to  the  nature  of  the  phase- change  front  and  can  be  extended  to 
complicated  multi- dimensional  problems  with  convective  terms  without 
involving  cumbersome  mathematical  schemes. 
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Y.  THERMAL  PROTECTION  FROI  INTENSE  LOCALIZED  IGVING  HEAT 
FLDIES  USING  PHASE- CHANGE  1ATERIALS 

5.1  Summary 

Various  technologies  which  protect  a  wall  from  being  burned  out  by 
an  intense  localized  moving  heat  source  have  been  reviewed,  and  a 
solution  to  this  problem  is  proposed  in  which  a  phase- change  material 
(PCM)  is  placed  underneath  the  wall  to  absorb  the  high  heat  flux.  The 
three-dimensional  melting  problem  is  nondimensionalized  and  modelled  with 
a  new  enthalpy  transforming  scheme.  The  numerical  results  show  that  the 
proposed  solution  reduces  the  peak  wall  temperature  significantly.  The 
method  of  coating  the  PCM  on  the  wall  surface  is  also  employed,  which  can 
maintain  the  surface  temperature  below  the  melting  temperature  of  the 
PCM. 
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5.2  Introduction 


Heat  Transfer  in  a  body  resulting  from  a  moving  heat  source  is  a 
well-known  heat  transfer  problem  and  has  had  broad  applications  in  the 
fields  of  welding,  laser  and  electron  beam  surface  processing  (Eckert  and 
Drake  [1972],  Festa  et  al.  [1988]),  etc.  In  these  cases,  a  surface  is 
subjected  to  a  localized  heat  input,  and  the  incoming  heat  flux  moves 
relatively  to  the  surface.  The  interest  for  the  above  applications  is 
how  to  heat  the  surface  effectively.  However,  in  some  applications  a 
surface  is  hit  by  a  moving  high  intensity  heat  source,  as  shown  in  Fig. 
5.1.  The  major  concern  here  is  how  to  protect  the  surface  from  being 
burned  out  by  the  moving  heat  flux.  This  is  indeed  a  problem  of  concern 
in  laser  thermal  threats  and  re-entry  situations  and  has  been  one  of  the 
motivations  of  the  present  study. 

The  methods  used  to  protect  surfaces  from  being  burned  out  by  a  high 
heat  flux  are  usually  ablation  and  heat  pipe  technologies.  Most  of  the 
previous  studies  on  ablation  concentrated  on  a  stationary  heat  input  with 
a  large  heating  area.  The  major  concern  for  this  problem  is  the  total 
energy  which  must  be  absorbed  during  a  given  time  period.  In  the 
ablation  technology  which  is  used  in  missile  re-entry  situations,  the 
body  surface  is  coated  with  a  solid  material  which  is  exposed  to  the  high 
heat  flux  and  is  allowed  to  melt  and  blow  away.  Thus,  a  large  amount  of 
the  incoming  heat  is  expended  in  melting  the  material  rather  than  being 
conducted  into  the  interior  of  the  vehicle,  so  that  the  temperature  of 
the  vehicle  surface  is  reduced.  The  thickness  of  the  coated  material 
must  be  greater  than  the  melted  thickness  during  this  time  period.  The 
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drawback  of  this  technology  is  that  once  the  coated  material  is  melted, 
it  is  immediately  blown  away.  The  coated  surface  can  be  used  only  once, 
and  for  the  next  mission  we  must  coat  the  surface  again.  For  vehicles 
which  have  frequent  missions  such  as  airplanes,  this  situation  may  be 
intolerable.  Moreover,  high  speed  vehicles  require  a  smooth  surface  to 
reduce  the  flow  resistance  and  to  reject  the  incoming  heat;  coated  or 
ablated  surfaces  may  have  difficulties  in  reaching  this  goal. 
Nevertheless,  due  to  the  simplicity  and  effectiveness  of  this  method,  it 
still  remains  an  alternative  to  protecting  the  surface  and  may  be  adopted 
to  the  case  of  a  moving  heat  source. 

Heat  pipe  technology  is  a  good  means  of  protecting  a  surface  from 
attack  by  a  high  heat  flux.  For  the  present  moving  heat  source  problem, 
however,  the  heat  pipe  may  only  work  under  melting  or  free  molecular 
conditions.  Furthermore,  a  heat  pipe  is  an  integral  vessel  which  cannot 
allow  vapor  leakage,  which  means  that  the  present  manufacturing 
technology  would  not  be  easily  adapted  to  vehicles  that  have  a  large 
surface  area.  For  these  reasons,  the  adoption  of  heat  pipe  technology 
for  the  present  problem  is  not  economical. 

Recently,  the  study  of  phase- change  materials  is  active  due  to 
applications  to  space- based  power  plants  and  the  utilization  of  solar 
energy.  Phase- change  materials  (PCM)  have  a  large  melting  heat,  so  it  is 
an  efficient  way  to  absorb  the  heat  energy  during  the  time  period  when 
the  materials  are  subject  to  heat  input  and  to  release  it  afterward  at  a 
relatively  constant  temperature.  It  is  advantageous  to  incorporate  the 
merits  of  the  above  technologies  and  propose  another  alternative  to 
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protect  the  surface  from  attack  by  a  high  heat  flux  for  the  present 
problem,  as  shown  in  Fig.  5.2.  There,  the  incoming  heat  input  moves 

along  the  surface  with  speed  U  and  the  heat  is  conducted  through  the 
outside  wall  to  the  PCM.  The  PCM  beneath  the  surface  melts  and  absorbs  a 
large  amount  of  the  incoming  heat.  Because  of  the  large  melting  latent 
heat  of  the  PCM  and  the  constant  melting  temperature  T  ,  the  peak  wall 
temperature  will  be  maintained  at  a  temperature  moderately  higher  than 
T  .  With  a  low  or  moderate  T  ,  the  reduction  of  the  peak  wall 
temperature  is  evident.  The  dividing  sheet,  the  soft  insulating 
material,  and  the  supporting  plate  may  be  used  to  prevent  the  separation 
of  the  PCM  from  the  surface  wall  during  the  melting  process,  and  to 
prevent  the  incoming  heat  from  being  conducted  into  the  cabin.  Because 
the  density  of  liquid  PCM  is  different  from  that  of  solid  PCM,  the  solid 
PCM  tends  to  separate  from  the  outside  wall,  or  to  form  voids  near  the 
wall,  which  reduces  the  heat  transfer  into  the  PCM.  If  the  thickness  of 
the  PCM  is  thin  and  the  density  difference  is  not  large,  this  situation 
will  not  be  very  severe. 

The  method  proposed  above  imposes  few  difficulties  on  the  present 
manufacturing  technology  and  might  also  be  used  to  cool  the  leading  edges 
of  space  vehicles  in  re-entry  situations.  Since  the  re-entry  time  is 
short,  a  moderate  thickness  of  the  PCM  will  greatly  reduce  the 
temperature  jump  at  the  leading  edges. 

Section  5.3  gives  a  nume’ical  analysis  of  this  problem  and  discusses 
the  important  parameters  to  protect  surfaces  from  an  intense  localized 
moving  heat  flux  using  PCM. 
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BY  HIGH  HEAT  FLUXES 


5.3  Analysis  of  Phase  Change  with  a  loving  Heat  Scarce 


The  present  problem  mentioned  previously  is  a  three-dimensional 
phase- change  problem  with  a  moving  heat  source.  Since  the  focus  here  is 
on  space  applications  and  the  PCM  thickness  is  comparatively  thin,  the 
natural  convection  due  to  gravity  will  not  be  considered.  Also,  since 
the  layer  of  the  PCM  is  thin,  the  influence  of  the  density  change  between 
the  solid  and  the  liquid  phases  is  expected  to  be  small,  and  will  be 
ignored. 

In  a  Cartesian  coordinate  system  (x',  y',  z')  fixed  at  the  outside 
wall  as  shown  in  Fig.  5.3,  the  heat  source  moves  relatively  to  the  plate 

with  speed  €  on  the  surface.  The  appropriate  energy  equation  is  (Kays 
and  Crawford  [1980]) 

d  (k £L)  _i(kOT)  d  (kOT 

dt  dx'  dx'  +  dy'  dy'  +  dz'  dz' 

where  the  relation  between  E  and  T  is  given  by  the  equation  of  state, 

HF 

_  =  c.  This  equation  is  applicable  to  the  three  different  regions, 

dT 

namely,  the  wall,  the  solid  PCM,  and  the  liquid  PCM  with  different 
relations  between  the  enthalpy  E,  the  temperature  T,  and  the  appropriate 
physical  properties. 

An  analysis  in  a  fixed  coordinate  system  is  difficult  for  this 
problem.  It  is  convenient  to  study  it  in  a  moving  coordinate  system 
where  the  origin  is  fixed  at  the  heated  spot.  Imagine  an  observer  riding 


(5.1) 
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along  with  the  incoming  heat  beam.  The  outside  wall  and  the  PCM  will 

-4 

travel  by  at  the  same  speed  -  C.  If  we  fix  a  moving  coordinate  system 
(x,  y,  z)  to  the  center  of  the  beam,  the  system  will  appear  in  reference 
to  the  fixed  coordinate  system  as  shown  in  Fig.  5.3. 


In  the  moving  coordinate  system,  the  problem  becomes  a  convective 
and  diffusive  heat  transfer  problem.  The  governing  equation  according  to 
Kays  and  Crawford  [1980]  is 

DE  d  5Tx  8  Ms  d  Ms 

P  _  =  _  (k  _ )  +  _  (k  _ )  +  (k 
Dt  dx  dx  dy  dy  dz  dz 


With  u  =  -U,  v  =  w  =  0,  we  have 


P 


DE  <9E 

_  =  /> 

Dt  dt 


p\} 


dl 

dx 


(5.3) 


Equation  (5.2)  then  becomes 

5E  0  3E  d  (k  M^  d  ^  0Tj  d  Mj 
M  di^di  St  fry  *  Jz  Tz 


(5.4) 


The  second  term  on  the  left-hand  side  of  the  equation  is  a  convective 
term,  while  the  terms  on  the  right-hand  side  are  diffusive  terms. 


The  relations  between  the  variables  of  the  two  coordinate  systems 
are  evident,  i.e.,  y  =  y,  z  -  z'  and  x  =  x  -  \\  U  dt.  When  U  is  a 
constant,  x  =  x'  -  Ut. 
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In  the  case  of  constant  physical  properties  for  the  wall,  the  solid 
PCM  and  the  liquid  PCM,  the  phase  change  occurs  at  a  single  temperature. 
Equation  (5.4)  can  be  written  for  each  region  as 


c 

w 

dT 

c 

w 

U^T  _ 

dn 

d2T 

dn 

r 

w 

F 

w 

dx 

w  + 

dy2  + 

Cs 

dT 

Cs 

dn 

dn 

dn 

F 

S 

dt 

fs 

dx 

^2  + 

w + 

$z2 

h 

dT 

C£ 

vdT 

d2T 

d2T 

dn 

V 

dt  ' 

5 

dx 

d x*  * 

W1  + 

dz2 

(wall) 


(solid  PCM)  (5.5) 


(liquid  PCM) 


where  C  =■  c  p,,  C  =  c  p  and  C,  =  c fpr  The  initial  and  boundary 

W  WW’S  oS  c  C 

conditions  are 


t  -  0,  T  -  T  ^ ,  0  <  z  <  h,  -od  <  x  <  cb,  -®<y<cD 


t  >  0 

dT 

k«  3i  =  V  2  =  °'  yI  -  Rh 


dT  ,  , _ 

kw  dx  =  *cT  ’  z  =  °’  +  y2  >  \ 


q  I  ^  =  q  |  ^  ,  z  =  6,  -  od  <  x  <  od  ,  ®  <  y  <  ®  (5-6) 


=0,  z  =  h,  -x  <  x  <  ®,  -x  <  y  <  m 
dz 
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T  =  T.,  0  <  z  <  h,  | x |  -*  od,  |y|  -*  ® 


u  =  U,  w  =  v  =  0,  0  <  z  <  h,  -®  <  x  <  ®,  -oD<y<oo 


At  the  liquid  and  solid  interface  within  the  PCM,  we  have 


Ts  =Tf 


St.  St 

1  w  +  s  w 


Hv  p 
nF 


(5.7) 


To  reduce  the  number  of  parameters  that  have  to  be  specified  for  the 
numerical  solutions,  the  following  dimensionless  variables  are 
introduced: 


* 


x 


* 

t 


C  R,  2 
w  h 


(5.8) 


Equations  (5.5- 5.7)  are  nondimensionalized  as 


<?T*  ffl*  d2T*  d2T*  <?2T* 

dt  dx  dx  ^  dy  ^  dz  ^ 


(wall) 


C  k  ry~* 
S  w  (IT 

C  k  * 
w  s  dt 


r  k  * 

S  Kw  L,+  5T 

C  k  T* 

w  s  ox 


<?2T*  d*T*  d2T* 

Vo  4  “  to  +  -  ~Vo 

f\  i. i  O  O 

ox  o  y  a  z 


(solid  PCM) 


(5.9) 
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(liquid  PCM) 


Cl  kw  W* 

C /  * 

1  w  u 

* 

fft 

i?2T* 

*  * 
52T  #2T 

C  k .  «*  C  k# 
w  l  at  w  l 

II 

■y 

* 

♦ 

* 

* 

* 

t  =  0, 

T  =  0,  0 

<  z  <  h  , 

-®  <  x  <  ®, 

* 

t  >  0 


^  ,  2*  =  0 


x*2  +  y*2  <  1 


m  1'  w 


Rt  (T  -  T.)3  at  *  *o  *2 

91  h  m  1  (T  +  1  )4,  z  =  0,  ,  x  2  +  y  L  >  1  (5.10) 

8z  \  m  i 


q|^*  =  q|^*+,  2  ~  ^  =  —  ,-®<x  <  ®,  -cd  <  y  <  ® 

Rh 


ffl  *  *  h  *  * 

- *■  =  0’  Z  =  h  —  S  J  -00<X  <01,  -CD  <  y  <  ® 

dz  h 


*  *  *  ,  *  I  1*1 

T  =  0,  0  <  z  <  h  ,  |x  |  -»  »,  |y  I  -*  b 

*  *  ♦  *  *  *  *  * 

u=U,  v  =  w  =  0,  0  <  z  <  h  ,  -qd<x<®,  -®  <  y  <  ® 


At  the  interface  of  the  liquid  and  solid  in  the  PCM 


*  * 

TS  ■  T<  ■  1 


cs  (Ta  •  Ti)  ,  k<  "t  ks  5TS 


(5.11) 


,  k<  "I  ks  ,  Cs  ah  v 

*v - W  +  j7~  W~  -  V  1 

v  dn  Kw  5n  w 
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The  independent  dimensionless  parameters  in  addition  to  those  in  eq. 
(5.8)  above  are 


Cs  ks  Vl  qh  Rh 

^  ^  ^  ’  <Tm  -  Ti)  kw 


*  *  R,  (T  -  T.)3 

c  u  h  v  m  1' 

)  *  »  h  , _ o  £  , 

k 


cs  '  Tr> 


T  *  T. 

m  l 


C  U  Rh  q.  Rh 

Upon  combining  w  with  !l  ’ 

k  (T  -  T-)  k~ 

w  '  m  w 


,  we  get  a  new  dimensionless 


parameter 


U  C  (T  -  T.) 
w  '  m  i' 


which  may  be  used  to  replace  the 


-T 

dimensionless  number  (^UR^/  k^.  Also,  since  h  -6  +  6  ,  the 

*  * 

independent  variable  h  may  be  replaced  by  b  .  Thus,  we  have 


T  "  f  (Nv>  Nq>  St,  1  >  Nr’  Nt’  Cs»-  CTu’  Ksv’  ^  >  ^p’  x  ’  1  ’  z  ) 


(5.12) 


where  T  =  (T  -  Tj)/^  -  T.) 

\  ■  UCw  Om  '  Tj )l\ 

\  -  %  V  <T«,  -  Ti>  k» 

St  =  cs  (T.  ■  Tj)/H 

t*  -  !<w  t/(CB  V) 

N  =  R,  (T  -  T.)3  a  el  k 
r  h  v  m  i'  '  w 

N.  =  T./(T  -  T.) 

t  i'  v  m  i; 

Cs»  -  Cs/Cw 

C(u  =  <VC» 

“sw  =  ks/k. 
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K£v  =  Vkv 

6*  --  6/ Rh 

6*  =  «n/R. 

^  P  h 

x  x/Rh 

y*  =  y/Rh 

z  z/Rh 
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5.4  Numerical  lodel  and  Solution  Procedure 


This  problem  is  nonlinear  due  to  the  moving  front,  so  that  exact 
analytical  solutions  for  this  type  of  nonlinear  problems  are  available 
only  for  some  simplified  and  idealized  systems.  Numerical  methods  appear 
to  be  the  only  practical  approach  for  handling  this  three-dimensional 
melting  problem  with  a  moving  heat  source. 

The  numerical  methods  used  to  solve  phase- change  problems  might  be 
divided  into  two  main  groups.  The  first  group  is  called  strong  numerical 
solutions;  examples  are  given  by  Okada  [1984],  and  Ho  and  Chen  [1986]. 
These  methods  are  only  applicable  to  processes  involving  one  or  two  space 
dimensions.  The  second  group  is  called  weak  numerical  solutions 
(Shamsunder  and  Sparrow  [1975],  Voller  and  Cross  [1981],  Cao  et  al. 
[1988],  Solomon  et  al.  [1986]).  These  methods  allow  us  to  avoid  paying 
explicit  attention  to  the  nature  of  the  phase- change  front.  Recently, 
Cao  et  al.  [1988]  proposed  a  new  enthalpy  transforming  model  which  starts 
from  eq.  (5.4)  and  transforms  the  equation  into  a  nonlinear  equation  with 
a  single  dependent  variable  E.  The  existing  algorithms  were  modified  to 
facilitate  the  numerical  calculations.  The  numerical  schemes  are 
flexible  and  can  handle  three-dimensional  problems  without  difficulty. 

To  simulate  the  circular  heat  source,  eq.  (5.4)  needs  to  be 
transformed  into  the  form  for  the  cylindrical  coordinate  system. 
Referring  to  Fig.  5.4  and  introducing  x  =  r  cos  0,  y  =  r  sin  0,  we  get 
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(b)  Top- View  of  the  Computational  Domain  of  the 
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(k  "  )  (5.13) 
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with  vr  =  -D  cos  0,  =  U  sin  0,  vz  =  0. 

Incorporating  the  continuity  equation 


1  8 
r  Hr 


(r  P  vr) 


1 

*  r  ~ST~  4  3t  =  0 


(5.14) 


we  obtain 


dpi  1  5(rVE)  1  5(V/E)  l  d  /  rkOT  \  1  d  ,  k5T  v 

dO  ~  t  dr  dr  +  r  JO  rdO  + 
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4  (kOT) 

"5“  “5"“ 

0Z  OZ 


(5.15) 


dE 

with  the  state  equation  _  =  c(T).  Following  the  analysis  given  by  Cao 

dT 

et  al.  [1988],  the  transforming  method  is  described  as  follows. 


In  the  case  of  constant  specific  heats  for  each  phase  and  that  the 
phase  change  occurs  at  a  single  temperature,  the  relation  between 
temperature  and  enthalpy  can  be  expressed  as  follows 


T  +  E/c 
m  '  s 


E  <  0 


0  <  E  <  H 


[Tm  +  (E  -  l)/ct  E  > 


(solid  phase) 

(mushy  phase)  (5.16) 

(liquid  phase) 
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Here,  E  =  0  has  been  selected  to  correspond  to  phase- change  materials  in 
their  solid  state  at  temperature  T  . 


The  "Kirchoff"  temperature  is  introduced  as 


T*  =  k(,)  d,  = 

m 


fk  (T  -  T  ) 
s  v  m; 

0 

■k<  (T  -  TJ 


1  < T. 
1R 

T  =  T 

m 

T  >  t 

m 


(5.17) 


Transforming  eq.  (5.16)  vith  the  definition  given  in  eq.  (5.17)  results 
in 


T+  - 


'  ks  E/cg  E  <  0 

0  0  <  E  <  H 

.  (E  -  H)/c^  E  >  H 


(5.18) 


An  enthalpy  function  is  introduced  as  follows 


T+  =  T(E)  E  +  S(E) 


(5.19) 


Eq.  (5.15)  is  transformed  into  a  nonlinear  equation  with  a  single 
dependent  variable  E. 
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5(E  p)  1  5(rVE^  1  5(V^E)  H  (r  5(rE)  ^  1  d  {  1  5(TE)  ^ 

+  r  5r  +  r  ffl  ~  t  Hr  5r"~  +  r  dd  r  50 


d  ,  9(rE)  , 

+  Tx  W  +  p 


(5.20) 


where 


1  d 


dS 


f  =  x  v  [r 

?  3F  3F 


1  5  rl  <9S  32S 

r  Jd  r  ffl  +  3z2 


v  =  -U  cos  9,  \g  =  U  sin  9 


r  =  r(E)  = 


and 


S  =  S(E)  = 


k  /c 
s'  s 

0 

-Vc^ 

0 

0 

* H  Vc£ 


E  <  0 

0  <  E  <  H 
E  >  H 

E  <  0 

0  <  E  <  H 
E  >  H 


For  the  wall,  no  phase  change  occurs  and  therefore 


(5.21) 


(5.22) 


E  -  c  T,  P  ?  k  /c  ,  and  S  e  0 
w  ’  w'  w’ 


(5.23) 
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with  eq.  (5.20)  still  being  applicable. 

The  initial  and  boundary  conditions  for  the  numerical  domain  in  Fig.  5.4 
are  as  follows: 


I 


s 


0  <  r  <  od,  0  <  0  <  360° 

0  <  r  <  uu,  0  <  0  <  360° 

r  <  Rh,  0  <  9  <  360° 

r  >  Rh,  0  <  9  <  360° 

0  <  r  <  0  <  0  <  3600 

0  <  r  <  ®,  0  <  0  <  360° 

r  =  RQ,  -90°  <  0  <  90« 

r  =  RQ,  90°  <  8  <  270° 

r  =  R  ,  90°  <  0  <  270° 


The  radius  RQ  must  be  sufficiently  large  such  that  the  region  r  > 
R05  '90°  <  0  <  90°  is  unaffected  by  the  moving  heat  source.  The  last  two 
boundary  conditions  imply  that  the  upwind  scheme  is  used,  and  the 
diffusive  term  is  neglected  for  the  outflow  boundary. 


Equation  (5.20)  is  discretized  with  the  numerical  scheme  proposed  by 
Cao  et  al.  [1988],  which  is  based  on  the  control- volume  finite-difference 
approach  described  by  Patankar  [1980].  The  general  three-dimensional 
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discretization  equation  for  the  phase- change  problem  with  a  fully 
implicit  scheme  in  the  cylindrical  coordinate  system  is 

apEp  =  agEg  +  ayEy  +  +  agE§  +  aTFT  +  aEFE  +  ^  (5.24) 


where 

AV 

aP  =  aPE  +  aPV  +  aPN  +  aPS  +  aPB  +  aPT  + 

aE  =  rEDe  +  max  f'Fe’  °]  >  aPE  =  rPDe  +  max  ^‘Fe’ 
ay  =  ryDw  +  max  [Fy,  0],  apy  =  rpDy  +  max  [Fw,  0] 

aN  =  FN^n  +  inaA  ^  Fn’  aPN  =  FP^n  +  max  ^  Fn’  ^ 
aS  =  TgDg  +  max  [Fg,  0],  apg  =  rpDg  +  max  [Fg,  0] 

aT  =  FTDt  +  max  t"Ft’  °J »  aPT  "  FPDt  +  max  t'Ft’ 

aB  =  rBDb  +  max  tFb’  °3»  aPB  =  FPDb  +  max  fFb’ 

Atyft 

b  =  It —  ^  +  De^E  Sp^  °w^Sp  Sy^  +  Dn^N  Sp^  Ds^P  ^  + 
Dt(Si  -  Sp)  -  D^(Sp  -  Sg) 


where  AV  is  the  volume  of  the  control  volume,  A 9  Az  (r  +  r  )/2,  the 
greater  of  a  and  b  is  given  by  max  [a,b],  E°  denotes  the  old  value  (at 
time  t)  of  E  at  grid  point  p,  and  denotes  the  old  value  of  p  at  grid 
point  p.  The  flow  rates  and  conductances  are  defined  as 


Fe  =  0>Ve  4r4z’  De 
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Because  of  the  nonlinearity  of  the  above  equation  and  the  implicit 
nature  of  the  scheme,  iterations  are  needed  at  each  time  step.  This 
procedure  is  the  same  as  those  which  solve  a  normal  nonlinear  equation. 


Since  the  dimensionless  parameters  in  the  last  section  are 
independent  of  the  coordinate  system,  eq.  (5.12)  is  still  applicable  to 
the  problem  in  the  cylindrical  coordinate  system  of  this  section. 
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5.5  Numerical  Results  and  Discussion 


The  discretization  equations  above  are  solved  with  the  Gauss- Seidel 
method.  The  grid  size  employed  is  32  x  30  x  20  with  the  grid  size  near 
the  wall  surface  being  finer.  To  check  the  validity  of  the  computer 
program,  the  calculation  has  been  made  with  a  moving  heat  source  without 
phase  change,  i.e.,  a  pure  wall  without  PCM,  and  the  steady- state  surface 
temperature  along  the  X-X  plane  (the  X-X  plane  corresponds  to  the  plane 
defined  by  $  =  0°  and  0  ~  180°  in  cylindrical  coordinates  or  y  =  0  in 
Cartesian  coordinates)  was  compared  with  the  analytical  result  of  a  point 
source  moving  on  the  surface  of  a  semi- inf inite  medium  by  Eckert  and 
Drake  [1972],  as  shown  in  Fig.  5.5.  The  steady- state  analytical  result 
has  the  form 

T  -  Ti  =  (Q/2xkwr)e'(U/2aw)  (r  + 

The  above  equation  was  obtained  from  a  moving  Cartesian  coordinate 
system  (x,  y,  z)  with  the  point  heat  source  at  the  origin  and  r  = 

v/x2  +  y2  +  z2. 

The  discrepancy  between  the  two  results  near  x  =  0  is  due  to  the 
nature  of  the  heat  source  that  was  modeled  in  each  case.  The  analytical 
result  has  an  infinitesimal  heat  source  at  x  =  0,  while  the  numerical  one 
has  a  heat  source  of  finite  radius  R^.  Considering  this,  the  accuracy  of 
the  numerical  results  is  very  good.  The  numerical  results  with  radiation 
(f  =  1)  into  space  for  the  surface  other  than  the  heated  spot  is  also 


(5.25) 
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FIG.5.5  COMPARISON  BETWEEN  ANALYTICAL  AND  NUMERICAL  RESULTS  WITHOUT  PHASE  CHANGE 


included  in  Fig.  5.5.  The  radiation  plays  no  significant  role  in  the 
temperature  distribution  of  the  wall  near  the  heat  source,  and  therefore 
will  not  be  included  in  the  next  discussion,  and  the  radiation  parameters 
Nf  and  will  be  dropped. 

Figure  5.6  shows  the  general  temperature  distribution  trend  on  the 

surface  as  a  three-dimensional  plot.  The  temperature  falls  off  sharply 

*  * 
in  the  portion  of  x  >0,  while  it  falls  off  gradually  for  x  <  0.  Also, 

the  symmetry  of  the  temperature  distribution  to  the  X-X  plane  indicated 

in  Fig.  5.4  is  evident.  Therefore,  the  representation  of  the  temperature 

distributions  will  focus  on  the  X-X  plane  for  simplicity. 

5.5.1  Moving  heat  source  without  phase  change 

As  a  first  step,  a  pure  wall  without  PCM  is  studied.  In  this 
situation,  eq.  (5.12)  reduces  to 

T  =  f(Nv,  Nq>  1-  »  8  ,  X  ,  y  ,  z  )  (5.26) 

Since  the  problem  is  now  independent  of  T  ,  (T  -  T.)  is  chosen  as  unity 

m  '  m  i'  J 

for  convenience. 

Figure  5.7  shows  the  variation  of  the  steady- state  wall  surface 
temperature  with  different  values  of  Ny.  A  larger  reduces  the  peak 
wall  temperature  significantly,  while  a  small  may  result  in  an 

intolerably  high  peak  wall  temperature.  The  dimensionless  number  Ng  has 
an  opposite  effect  on  the  wall  temperature,  as  shown  in  Fig.  5.8.  A 
smaller  Ng  corresponds  to  a  lower  peak  wall  temperature,  while  a  larger 
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FIG. 5.8  GRAPH  OF  THE  SURFACE  TEMPERATURE  DISTRIBUTION  TREND 
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X*  ALONG  THE  X-X  PLANE  AT  1=0 

FIG.5.8  STEADY  STATE  TEMPERATURE  PROFILES  WITH  DIFFERENT  NQ  WITHOUT  PCM 


Nq  corresponds  to  a  higher  peak  wall  temperature.  Ve  can  conclude  that 

* 

with  a  certain  6  ,  when  is  large  enough  or  Nq  is  small  enough,  the 

peak  wall  temperature  is  low  and  the  wall  does  not  need  to  be  protected. 

When  Ny  is  small  or  Nq  is  large,  something  must  be  done  to  protect  the 

* 

wall,  otherwise  the  wall  will  be  burned  out.  Also,  a  larger  6  has  the 
effect  of  relieving  the  peak  wall  temperature,  but  is  often  not  practical 
in  space  applications. 

5.5.2  Moving  heat  source  with  PCM  underneath  the  wall 

As  mentioned  before,  one  of  the  alternatives  of  protecting  the  wall 

is  to  put  the  PCM  beneath  the  wall  as  shown  in  Figs.  5.2  and  5.4.  Figure 

* 

5.9  shows  the  numerical  isotherms  of  the  dimensionless  temperature  T  = 

(T  -  T^)/  (Tm  -  'T)  for  t  =  0.1  second  along  the  XX  plane.  Other 

parameters  are  Ny  =  1.2,  Nq  =  60,  St  =  0.001,  Cgw  =  0.94,  C ^  =  1.38,  K ^ 

*  * 

=  K  =  0.65,  6  -  0.51,  and  h  =  1.45.  The  center  of  the  heat  source  is 
located  at  x  =  0.0.  The  solid  line  labelled  with  1.0  indicates  the 
melting  front  at  this  time.  After  about  0.5  second,  the  steady- state 
condition  is  reached.  Figure  5.10  shows  the  corresponding  steady- state 
isotherms  of  the  dimensionless  temperature  at  the  same  plane.  The 
melting  front  is  also  labelled  with  1.0. 

Figure  5.11  shows  the  steady  state  wall  surface  temperature  along 
the  X-X  plane  with  different  Stefan  numbers,  St  -  cg(Tm  -  T^)/H.  Also 
shown  in  the  figure  is  that  of  a  pure  wall  (without  PCM)  having  the  same 

weight  per  unit  surface  area  as  that  of  the  wall- PCM  module.  The 

reduction  of  the  peak  wall  temperature  is  significant.  In  the  figure, 

* 

^m  =  ^m^h  *s  dimensionless  maximum  melting  front  depth,  where  8m  is 
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FIG.5.9  ISOTHERMS  OF  THE  SOLUTION  for  t-O.ls  AT  X-X  PLANE 
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X*  ALONG  THE  X-X  PLANE  AT  2-0 

.11  STEADY  STATE  TEMPERATURE  PROFILES  WITH  DIFFERENT  St 


the  maximum  melting  front  depth  from  the  wall  surface  corresponding  to 
the  steady  state. 


* 

The  dimensionless  wall  thickness  6  has  a  significant  influence  on 

* 

the  temperature  distribution.  With  a  smaller  6  the  reduction  of  the 
peak  wall  temperature  is  more  evident,  as  shown  in  Fig.  5.12. 

Like  the  moving  heat  source  problem  without  phase  change,  when  Ny  is 
small  enough  and  Nq  is  large  enough  the  dimensionless  temperature  will 
become  high,  as  shown  in  Fig.  5.13.  The  case  with  Ny  =  0.3  and  Ng  =  60 
has  a  peak  wall  temperature  almost  four  times  as  high  as  that  with  Ny  = 
2.4  and  Ng  =  30.  In  this  situation,  perhaps  we  have  to  resort  to  another 
alternative. 

5.5.3  Moving  heat  source  with  TZ'A  coated  on  the  surface 

The  configuration  of  this  technology  of  protecting  the  surface  is 
illustrated  in  Fig.  5.14.  The  numerical  procedure  is  similar  to  that 
with  PCM  beneath  the  wall.  The  calculation  is  conducted  with  three 
different  Stefan  numbers  and  the  results  are  shown  in  Fig.  5.15.  The 
liquid  phase  in  the  numerical  domain  is  assumed  not  to  be  blown  away  in 
this  situation.  Also  shown  in  the  tigure  is  the  temperature  distribution 
of  a  pure  plate  (without  PCM)  having  the  same  weight  per  unit  surface 
area  as  that  of  the  PCM- wall  module  indicated  in  Fig.  5.14.  It  is 
evident  that  the  dimensionless  temperature  of  the  wall  surface  with  PCM 
coated  on  it  has  no  way  to  exceed  one,  which  means  that  the  wall  surface 
temperature  will  be  less  than  T  ,  provided  than  the  maximum  melting  front 
depth  is  less  that  the  thickness  of  the  PCM  coat  on  the  surface. 
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(a)  Side-View  of  the  Computational  Domain  of  the 

PCM  -  WALL  Module. 


(b)  Top- View  of  the  Computational  Domain  of  the 

PCM  -  WALL  Module. 

FIG.  5.H  PICTORIAL  DESCRIPTION  OF  THE  COMPUTATIONAL  DOMAIN  FOR 

THE  PCM  -  WALL  MODULE. 
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FIG. 5. 15  STEADY  STATE  TEMPERATURE  VARIATION  OF  WALL  SURFACE  WITH  St 


5.6  Conclusions  and  Remarks 


If  a  wall  surface  is  subject  to  an  intense  moving  heat  flux,  and 

with  a  large  Nq  and  a  small  Ny,  the  surface  is  prone  to  be  burned  out 

unless  some  measure  is  taken  to  protect  it.  The  use  of  PCM  proves  to  be 

a  good  means  of  protecting  the  surface  in  this  respect.  One  way  to 

protect  the  surface  is  to  put  the  PCM  underneath  the  wall.  This  method 

* 

is  more  efficient  with  smaller  St  and  S  .  With  increasing  and  Nv 
decreasing  further,  the  method  of  coating  the  PCM  on  the  surface  is 
needed  to  protect  the  surface.  In  this  situation,  the  wall  surface 
temperature  will  be  less  than  Tm  provided  that  the  maximum  melting  depth 
is  less  than  the  thickness  of  the  PCM  coat  on  the  surface.  Ve  should 
point  out  that  both  methods  have  advantages  and  disadvantages.  Even  if 
the  method  of  putting  the  PCM  underneath  the  wall  is  less  efficient  when 
Nq  and  Ny  are  beyond  some  limits,  it  is  ready  to  be  used  again  and  ^uts 
no  restrictions  on  the  wall  surface.  In  these  cases,  a  trade-off  will  be 
reached  when  selecting  the  protective  technology. 
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